{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "GREMLIN_TF_v2.ipynb",
      "version": "0.3.2",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/sokrypton/GREMLIN_CPP/blob/master/GREMLIN_TF.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "metadata": {
        "id": "mu2S86VhS-8t",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# GREMLIN_TF v2.1\n",
        "GREMLIN implemented in tensorflow\n",
        "\n",
        "### Change log:\n",
        "*   02Apr2019\n",
        " - fixing a few hard-coded values, to allow GREMLIN to work with any alphabet (binary, protein, rna etc)\n",
        "*   22Jan2019\n",
        " - moving [GREMLIN_TF_simple](https://colab.research.google.com/github/sokrypton/GREMLIN_CPP/blob/master/GREMLIN_TF_simple.ipynb) to a seperate notebook\n",
        "*   19Jan2019\n",
        " - in the past we found that optimizing V first, required less iterations for convergence. Since V can be computed exactly (assuming no W), we replace this first optimization step with a simple V initialization.\n",
        " - a few variables were renamed to be consistent with the c++ version\n",
        "*   16Jan2019\n",
        " - updating how indices are handled (for easier/cleaner parsing)\n",
        " - minor speed up in how we symmetrize and zero the diagional of W\n",
        "*   15Jan2019\n",
        " - LBFGS optimizer replaced with a modified version of the ADAM optimizer\n",
        " - Added option for stochastic gradient descent (via batch_size)\n",
        "  \n",
        "### Method:\n",
        "GREMLIN takes a multiple sequence alignment (MSA) and returns a Markov Random Field (MRF). The MRF consists of a one-body term (V) that encodes conservation, and a two-body term (W) that encodes co-evolution.\n",
        "\n",
        "For more details about the method see:\n",
        "[Google slides](https://docs.google.com/presentation/d/1aooxoksosSv7CWs9-ktqhUjyXR3wrgbG5a6PCr92od4/) and accompanying [Google colab](https://colab.research.google.com/drive/17RJcExuyifnd7ShTcsZGh6mBpWq0-s60)\n",
        "\n",
        "See [GREMLIN_TF_simple](https://colab.research.google.com/github/sokrypton/GREMLIN_CPP/blob/master/GREMLIN_TF_simple.ipynb) for a stripped down version of this code (with no funky gap removal, sequence weight, etc). This is intented for educational purpose,  and could also be very useful for anyone trying to modify or improve the algorithm!\n"
      ]
    },
    {
      "metadata": {
        "id": "yM3wyYU5SwYn",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# ------------------------------------------------------------\n",
        "# \"THE BEERWARE LICENSE\" (Revision 42):\n",
        "# <so@g.harvard.edu> and <pkk382@g.harvard.edu> wrote this code.\n",
        "# As long as you retain this notice, you can do whatever you want\n",
        "# with this stuff. If we meet someday, and you think this stuff\n",
        "# is worth it, you can buy us a beer in return.\n",
        "# --Sergey Ovchinnikov and Peter Koo\n",
        "# ------------------------------------------------------------"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "NLUvPVyxb7bo",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## libraries"
      ]
    },
    {
      "metadata": {
        "id": "YyJpLM_tJfrY",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# IMPORTANT, only tested using PYTHON 3!\n",
        "import numpy as np\n",
        "import tensorflow as tf\n",
        "import matplotlib.pylab as plt\n",
        "from scipy import stats\n",
        "from scipy.spatial.distance import pdist,squareform\n",
        "import pandas as pd"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "j0Yp7bPRmvwU",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## Params"
      ]
    },
    {
      "metadata": {
        "id": "o3c7KURqmugY",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "################\n",
        "# note: if you are modifying the alphabet\n",
        "# make sure last character is \"-\" (gap)\n",
        "################\n",
        "alphabet = \"ARNDCQEGHILKMFPSTWYV-\"\n",
        "states = len(alphabet)\n",
        "a2n = {}\n",
        "for a,n in zip(alphabet,range(states)):\n",
        "  a2n[a] = n\n",
        "################\n",
        "\n",
        "def aa2num(aa):\n",
        "  '''convert aa into num'''\n",
        "  if aa in a2n: return a2n[aa]\n",
        "  else: return a2n['-']"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "bX6GXKV3I2pm",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## Functions for prepping the MSA (Multiple sequence alignment)"
      ]
    },
    {
      "metadata": {
        "id": "zA0Bne59SUIu",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# from fasta\n",
        "def parse_fasta(filename,limit=-1):\n",
        "  '''function to parse fasta'''\n",
        "  header = []\n",
        "  sequence = []\n",
        "  lines = open(filename, \"r\")\n",
        "  for line in lines:\n",
        "    line = line.rstrip()\n",
        "    if line[0] == \">\":\n",
        "      if len(header) == limit:\n",
        "        break\n",
        "      header.append(line[1:])\n",
        "      sequence.append([])\n",
        "    else:\n",
        "      sequence[-1].append(line)\n",
        "  lines.close()\n",
        "  sequence = [''.join(seq) for seq in sequence]\n",
        "  return np.array(header), np.array(sequence)\n",
        "\n",
        "def filt_gaps(msa,gap_cutoff=0.5):\n",
        "  '''filters alignment to remove gappy positions'''\n",
        "  tmp = (msa == states-1).astype(np.float)\n",
        "  non_gaps = np.where(np.sum(tmp.T,-1).T/msa.shape[0] < gap_cutoff)[0]\n",
        "  return msa[:,non_gaps],non_gaps\n",
        "\n",
        "def get_eff(msa,eff_cutoff=0.8):\n",
        "  '''compute effective weight for each sequence'''\n",
        "  ncol = msa.shape[1]\n",
        "  \n",
        "  # pairwise identity\n",
        "  msa_sm = 1.0 - squareform(pdist(msa,\"hamming\"))\n",
        "\n",
        "  # weight for each sequence\n",
        "  msa_w = (msa_sm >= eff_cutoff).astype(np.float)\n",
        "  msa_w = 1/np.sum(msa_w,-1)\n",
        "  \n",
        "  return msa_w\n",
        "\n",
        "def mk_msa(seqs):\n",
        "  '''converts list of sequences to msa'''\n",
        "  \n",
        "  msa_ori = []\n",
        "  for seq in seqs:\n",
        "    msa_ori.append([aa2num(aa) for aa in seq])\n",
        "  msa_ori = np.array(msa_ori)\n",
        "  \n",
        "  # remove positions with more than > 50% gaps\n",
        "  msa, v_idx = filt_gaps(msa_ori,0.5)\n",
        "  \n",
        "  # compute effective weight for each sequence\n",
        "  msa_weights = get_eff(msa,0.8)\n",
        "\n",
        "  # compute effective number of sequences\n",
        "  ncol = msa.shape[1] # length of sequence\n",
        "  w_idx = v_idx[np.stack(np.triu_indices(ncol,1),-1)]\n",
        "  \n",
        "  return {\"msa_ori\":msa_ori,\n",
        "          \"msa\":msa,\n",
        "          \"weights\":msa_weights,\n",
        "          \"neff\":np.sum(msa_weights),\n",
        "          \"v_idx\":v_idx,\n",
        "          \"w_idx\":w_idx,\n",
        "          \"nrow\":msa.shape[0],\n",
        "          \"ncol\":ncol,\n",
        "          \"ncol_ori\":msa_ori.shape[1]}\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "fky6gk-HlFyi",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## GREMLIN"
      ]
    },
    {
      "metadata": {
        "id": "gx14M7Tvu-Ct",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# external functions\n",
        "\n",
        "def sym_w(w):\n",
        "  '''symmetrize input matrix of shape (x,y,x,y)'''\n",
        "  x = w.shape[0]\n",
        "  w = w * np.reshape(1-np.eye(x),(x,1,x,1))\n",
        "  w = w + tf.transpose(w,[2,3,0,1])\n",
        "  return w\n",
        "\n",
        "def opt_adam(loss, name, var_list=None, lr=1.0, b1=0.9, b2=0.999, b_fix=False):\n",
        "  # adam optimizer\n",
        "  # Note: this is a modified version of adam optimizer. More specifically, we replace \"vt\"\n",
        "  # with sum(g*g) instead of (g*g). Furthmore, we find that disabling the bias correction\n",
        "  # (b_fix=False) speeds up convergence for our case.\n",
        "  \n",
        "  if var_list is None: var_list = tf.trainable_variables() \n",
        "  gradients = tf.gradients(loss,var_list)\n",
        "  if b_fix: t = tf.Variable(0.0,\"t\")\n",
        "  opt = []\n",
        "  for n,(x,g) in enumerate(zip(var_list,gradients)):\n",
        "    if g is not None:\n",
        "      ini = dict(initializer=tf.zeros_initializer,trainable=False)\n",
        "      mt = tf.get_variable(name+\"_mt_\"+str(n),shape=list(x.shape), **ini)\n",
        "      vt = tf.get_variable(name+\"_vt_\"+str(n),shape=[], **ini)\n",
        "      \n",
        "      mt_tmp = b1*mt+(1-b1)*g\n",
        "      vt_tmp = b2*vt+(1-b2)*tf.reduce_sum(tf.square(g))\n",
        "      lr_tmp = lr/(tf.sqrt(vt_tmp) + 1e-8)\n",
        "\n",
        "      if b_fix: lr_tmp = lr_tmp * tf.sqrt(1-tf.pow(b2,t))/(1-tf.pow(b1,t))\n",
        "\n",
        "      opt.append(x.assign_add(-lr_tmp * mt_tmp))\n",
        "      opt.append(vt.assign(vt_tmp))\n",
        "      opt.append(mt.assign(mt_tmp))\n",
        "        \n",
        "  if b_fix: opt.append(t.assign_add(1.0))\n",
        "  return(tf.group(opt))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "BqYqlJAXVI9N",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def GREMLIN(msa, opt_type=\"adam\", opt_iter=100, opt_rate=1.0, batch_size=None):\n",
        "  \n",
        "  ##############################################################\n",
        "  # SETUP COMPUTE GRAPH\n",
        "  ##############################################################\n",
        "  # kill any existing tensorflow graph\n",
        "  tf.reset_default_graph()\n",
        "  \n",
        "  ncol = msa[\"ncol\"] # length of sequence\n",
        "\n",
        "  # msa (multiple sequence alignment) \n",
        "  MSA = tf.placeholder(tf.int32,shape=(None,ncol),name=\"msa\")\n",
        "  \n",
        "  # one-hot encode msa\n",
        "  OH_MSA = tf.one_hot(MSA,states)\n",
        "  \n",
        "  # msa weights\n",
        "  MSA_weights = tf.placeholder(tf.float32, shape=(None,), name=\"msa_weights\")\n",
        "  \n",
        "  # 1-body-term of the MRF\n",
        "  V = tf.get_variable(name=\"V\", \n",
        "                      shape=[ncol,states],\n",
        "                      initializer=tf.zeros_initializer)\n",
        "\n",
        "  # 2-body-term of the MRF\n",
        "  W = tf.get_variable(name=\"W\",\n",
        "                      shape=[ncol,states,ncol,states],\n",
        "                      initializer=tf.zeros_initializer)\n",
        "  \n",
        "  # symmetrize W\n",
        "  W = sym_w(W)\n",
        "  \n",
        "  def L2(x): return tf.reduce_sum(tf.square(x))\n",
        "  \n",
        "  ########################################\n",
        "  # V + W\n",
        "  ########################################\n",
        "  VW = V + tf.tensordot(OH_MSA,W,2)\n",
        "\n",
        "  # hamiltonian\n",
        "  H = tf.reduce_sum(tf.multiply(OH_MSA,VW),axis=(1,2))\n",
        "  # local Z (parition function)\n",
        "  Z = tf.reduce_sum(tf.reduce_logsumexp(VW,axis=2),axis=1)\n",
        "\n",
        "  # Psuedo-Log-Likelihood\n",
        "  PLL = H - Z\n",
        "\n",
        "  # Regularization\n",
        "  L2_V = 0.01 * L2(V)\n",
        "  L2_W = 0.01 * L2(W) * 0.5 * (ncol-1) * (states-1)\n",
        "\n",
        "  # loss function to minimize\n",
        "  loss = -tf.reduce_sum(PLL*MSA_weights)/tf.reduce_sum(MSA_weights)\n",
        "  loss = loss + (L2_V + L2_W)/msa[\"neff\"]\n",
        "\n",
        "  ##############################################################\n",
        "  # MINIMIZE LOSS FUNCTION\n",
        "  ##############################################################\n",
        "  if opt_type == \"adam\":  \n",
        "    opt = opt_adam(loss,\"adam\",lr=opt_rate)\n",
        "\n",
        "  # generate input/feed\n",
        "  def feed(feed_all=False):\n",
        "    if batch_size is None or feed_all:\n",
        "      return {MSA:msa[\"msa\"], MSA_weights:msa[\"weights\"]}\n",
        "    else:\n",
        "      idx = np.random.randint(0,msa[\"nrow\"],size=batch_size)\n",
        "      return {MSA:msa[\"msa\"][idx], MSA_weights:msa[\"weights\"][idx]}\n",
        "  \n",
        "  # optimize!\n",
        "  with tf.Session() as sess:\n",
        "    # initialize variables V and W\n",
        "    sess.run(tf.global_variables_initializer())\n",
        "    \n",
        "    # initialize V\n",
        "    msa_cat = tf.keras.utils.to_categorical(msa[\"msa\"],states)\n",
        "    pseudo_count = 0.01 * np.log(msa[\"neff\"])\n",
        "    V_ini = np.log(np.sum(msa_cat.T * msa[\"weights\"],-1).T + pseudo_count)\n",
        "    V_ini = V_ini - np.mean(V_ini,-1,keepdims=True)\n",
        "    sess.run(V.assign(V_ini))\n",
        "    \n",
        "    # compute loss across all data\n",
        "    get_loss = lambda: round(sess.run(loss,feed(feed_all=True)) * msa[\"neff\"],2)\n",
        "    print(\"starting\",get_loss())\n",
        "    \n",
        "    if opt_type == \"lbfgs\":\n",
        "      lbfgs = tf.contrib.opt.ScipyOptimizerInterface\n",
        "      opt = lbfgs(loss,method=\"L-BFGS-B\",options={'maxiter': opt_iter})\n",
        "      opt.minimize(sess,feed(feed_all=True))\n",
        "      \n",
        "    if opt_type == \"adam\":\n",
        "      for i in range(opt_iter):\n",
        "        sess.run(opt,feed())  \n",
        "        if (i+1) % int(opt_iter/10) == 0:\n",
        "          print(\"iter\",(i+1),get_loss())\n",
        "    \n",
        "    # save the V and W parameters of the MRF\n",
        "    V_ = sess.run(V)\n",
        "    W_ = sess.run(W)\n",
        "    \n",
        "  # only return upper-right triangle of matrix (since it's symmetric)\n",
        "  tri = np.triu_indices(ncol,1)\n",
        "  W_ = W_[tri[0],:,tri[1],:]\n",
        "  \n",
        "  mrf = {\"v\": V_,\n",
        "         \"w\": W_,\n",
        "         \"v_idx\": msa[\"v_idx\"],\n",
        "         \"w_idx\": msa[\"w_idx\"]}\n",
        "  \n",
        "  return mrf"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "mppg0JLtP25z",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## EXAMPLE"
      ]
    },
    {
      "metadata": {
        "id": "bnigmLmAlyWv",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# download example fasta MSA\n",
        "!wget -q -nc https://gremlin2.bakerlab.org/db/PDB_EXP/fasta/4FAZA.fas"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "osaZwTSMOicF",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# ===============================================================================\n",
        "# PREP MSA\n",
        "# ===============================================================================\n",
        "# parse fasta\n",
        "names, seqs = parse_fasta(\"4FAZA.fas\")\n",
        "\n",
        "# process input sequences\n",
        "msa = mk_msa(seqs)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "CoBRuqmVVrbD",
        "colab_type": "code",
        "outputId": "465a2aa2-8d31-4835-dc89-2103262262c5",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 360
        }
      },
      "cell_type": "code",
      "source": [
        "%%time\n",
        "# ===============================================================================\n",
        "# RUN GREMLIN\n",
        "# ===============================================================================\n",
        "# Note: the original GREMLIN uses the \"lbfgs\" optimizer which is EXTREMELY slow \n",
        "# in tensorflow. The modified adam optimizer is much faster, but may \n",
        "# require adjusting number of iterations (opt_iter) to converge to the same \n",
        "# solution. To switch back to the original, set opt_type=\"lbfgs\".\n",
        "# ===============================================================================\n",
        "mrf = GREMLIN(msa)"
      ],
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "WARNING:tensorflow:From /usr/local/lib/python3.6/dist-packages/tensorflow/python/framework/op_def_library.py:263: colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.\n",
            "Instructions for updating:\n",
            "Colocations handled automatically by placer.\n",
            "WARNING:tensorflow:From /usr/local/lib/python3.6/dist-packages/tensorflow/python/ops/math_ops.py:3066: to_int32 (from tensorflow.python.ops.math_ops) is deprecated and will be removed in a future version.\n",
            "Instructions for updating:\n",
            "Use tf.cast instead.\n",
            "starting 50463.14\n",
            "iter 10 29873.72\n",
            "iter 20 27907.74\n",
            "iter 30 26972.23\n",
            "iter 40 26687.46\n",
            "iter 50 26576.22\n",
            "iter 60 26540.93\n",
            "iter 70 26524.79\n",
            "iter 80 26518.27\n",
            "iter 90 26515.09\n",
            "iter 100 26513.35\n",
            "CPU times: user 1.68 s, sys: 521 ms, total: 2.2 s\n",
            "Wall time: 2.19 s\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "id": "jCrfC2Um4xww",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## Explore the contact map\n",
        "### Contact prediction:\n",
        "\n",
        "For contact prediction, the W matrix is reduced from LxLx21x21 to LxL matrix (by taking the L2norm for each of the 20x20). In the code below, you can access this as mtx[\"raw\"]. Further correction (average product correction) is then performed to the mtx[\"raw\"] to remove the effects of entropy, mtx[\"apc\"]. The relative ranking of mtx[\"apc\"] is used to assess importance. When there are enough effective sequences (>1000), we find that the top 1.0L contacts are ~90% accurate! When the number of effective sequences is lower, NN can help clean noise and fill in missing contacts.\n"
      ]
    },
    {
      "metadata": {
        "id": "XxgQArVUPyPH",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## Functions for extracting contacts from MRF"
      ]
    },
    {
      "metadata": {
        "id": "nMxp7up_P1_q",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "###################\n",
        "def normalize(x):\n",
        "  x = stats.boxcox(x - np.amin(x) + 1.0)[0]\n",
        "  x_mean = np.mean(x)\n",
        "  x_std = np.std(x)\n",
        "  return((x-x_mean)/x_std)\n",
        "\n",
        "def get_mtx(mrf):\n",
        "  '''get mtx given mrf'''\n",
        "  \n",
        "  # l2norm of 20x20 matrices (note: we ignore gaps)\n",
        "  raw = np.sqrt(np.sum(np.square(mrf[\"w\"][:,:-1,:-1]),(1,2)))\n",
        "  raw_sq = squareform(raw)\n",
        "\n",
        "  # apc (average product correction)\n",
        "  ap_sq = np.sum(raw_sq,0,keepdims=True)*np.sum(raw_sq,1,keepdims=True)/np.sum(raw_sq)\n",
        "  apc = squareform(raw_sq - ap_sq, checks=False)\n",
        "\n",
        "  mtx = {\"i\": mrf[\"w_idx\"][:,0],\n",
        "         \"j\": mrf[\"w_idx\"][:,1],\n",
        "         \"raw\": raw,\n",
        "         \"apc\": apc,\n",
        "         \"zscore\": normalize(apc)}\n",
        "  return mtx\n",
        "\n",
        "def plot_mtx(mtx,key=\"zscore\",vmin=1,vmax=3):\n",
        "  '''plot the mtx'''\n",
        "  plt.figure(figsize=(5,5))\n",
        "  plt.imshow(squareform(mtx[key]), cmap='Blues', interpolation='none', vmin=vmin, vmax=vmax)\n",
        "  plt.grid(False)\n",
        "  plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "RSleviAVPJ36",
        "colab_type": "code",
        "outputId": "b9842b7f-6878-460c-c34e-4a239f0aa9fe",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 321
        }
      },
      "cell_type": "code",
      "source": [
        "mtx = get_mtx(mrf)  \n",
        "plot_mtx(mtx)"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAATMAAAEwCAYAAADbzJbwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt8VNW5N/BfSBhDLkAImYRwEYoQ\nIheBwlu5RWMoNdHzEumHEqfAay1UpCC2L5eIqD2VA3IxitgeLAK+RY9EgSJVKlRaXkWHVKgiQQ6Y\nVJBLyAWHkJAZSIY5f/ApZe01sPfszEx21vy+f/ks1957gcnjzLPXJcrn8/lARNTKtWnpARARBQOT\nGREpgcmMiJTAZEZESmAyIyIlMJkRkRJizF64ZMkSHDx4EFFRUVi4cCEGDRoUzHEREQXEVDL729/+\nhhMnTqC4uBjl5eVYuHAhiouLb9jf06R/T9fFy0KcFG8zMzS6TpP3ihDHRMsfxHccrhDivP5dQjqm\nlnL2vEeI0zrG6l6jnYIZFRUV1DGRMbEGs5Spr5lOpxNjx44FAPTu3Ru1tbWor683cysioqAwlcxq\namqQlJR0Le7UqROqq6uDNigiokAF5QUAV0QRUUszVTOz2+2oqam5FldVVSElJaVZA9HWyP7+tUvq\nM7RXkhAfOX1B6pPZtb0QX7kiJ9o2bfRrHxc1hb54o1/cLcRfjUzLTI2spu6SEHdOvCXgexh1uUms\n+9lizP3/197+5mP869EqqS07wy7Egxa+L/X5Ysm9psYTCpN/f0Bqe33qd4X4RE2D1OfWznFBef7B\nE+eF+I5bOwblvkaZ+skYNWoUdu7cCQA4fPgw7HY7EhISgjowIqJAmPq4MXToUPTv3x8FBQWIiorC\nM888E+xxEREFxPR3p7lz5wZzHEREzcIVAESkhKhwbM5YURv4hNhT37qltg/KK4X4oeE9mzWu1uKr\ns/Icvj5pgdcoH3xtv9T25kPDbnpNbUOj1NZz3JNC7Nq7POCxAOYmslLkCemkWSIiq2EyIyIlMJkR\nkRLCUjPTLjT3NyHW3kGsl3Tr1E73vkk/XCO1ubbMEOIPj8nLrLL6Nm+CrxWd1yzUB4COXKxv2PHq\ni1Jbz5T4gO9z5+LdUtu+RTkBP9/Ms82yeu2SNTMiiihMZkSkBCYzIlICkxkRKaFFXgD489qnx4U4\nUibEGvH0+0eltl/fmxHwfZLuf0Fqc737i4Dv477sFeJ2tuiA79EavHPotBCPH9hV95qqWo/Upn25\nRYHhCwAiiihMZkSkBCYzIlJCi9TMjOwQa2RCrD+lJ2uFeED3DgZGSKrTO2mp3k9hN0FTrPns+Hmp\nz5Ce4m6q9im/l/pUbZwqxB+X1Uh9Rt3WWWqjq1gzI6KIwmRGREpgMiMiJTCZEZESLDNp1gxtsR+Q\nC/6NmqPKAKCt5rgy7RFZALDqk+NCvP7BwSZGKKu6IB7TpncEmhVpC+HaIrhRq/eWC/Hs0b1Nj6ml\nPL7tsBC/mN9f6rP3K7HgP7pPaIr9Fy/Jv2jxt7S+IxK1+AKAiCIKkxkRKYHJjIiU0CI1sytX5Ee2\naSNOYjS7Q6y2RqatjwHA9E0HhXhtwR1Sn6Nn6oQ4Iz1R99lkXU1e8eciJjrw/48bqb8Gy56j8s//\n3RnN3yH5jEs+9Sw9SdzV+ZuaBqlPj85xuvf+/qq9QvznOaMDHJ1/rJkRUURhMiMiJTCZEZESmMyI\nSAmtetKsEdpiPyAX/LdpdhQFgE63iMe0qXg8XWug3bnV7K6twXgBEEpJ9z0vxHt++1Opzx23mpuc\n3NrxBQARRRQmMyJSApMZESlBuZqZdtG4vzqDtkaW7+fUnYff/FyIg7XQnIgCw5oZEUUUJjMiUgKT\nGREpwVAyO3bsGMaOHYvXX38dAFBRUYEpU6bA4XBgzpw5uHz5ckgHSUSkR7e01tDQgGeffRYjRoy4\n1vbSSy/B4XAgNzcXRUVF2Lx5MxwOh+GHXvTzRiDeaJVPh3aH2AVt5d1LtRNitcV+gAV/q9DusKLd\nXcUf7QRZwHqTZCn4dP8L22w2rF27Fna7/VpbSUkJcnJyAADZ2dlwOp2hGyERkQG6H4diYmIQEyN2\nc7vdsNmufrpJTk5GdbW89xIRUTg1+7N3GKapERHpMlWoiouLg8fjQWxsLCorK4WvoEYEqz7mj5Fa\nVwbEXWODtYh82K/+LMTln3wq9VmzbLIQTxrSQ+pTUyee4NQ5sfWd4BQsRmpkWqyPRSZT/9VHjhyJ\nnTt3AgB27dqFMWPGBHVQRESB0v2IVFpaimXLluH06dOIiYnBzp07sXLlShQWFqK4uBjp6enIz88P\nx1iJiG5IN5kNGDAAGzdulNo3bNgQkgEREZmh3ELzllTb0CjEHeLa6l7z003yHLd1BaGZ49bg58Tr\nOJ0Tr89flCdEd4y3+elJoXbn4t1CvP0x8fQje/vQ1Va9mvl+8/54ROpTNP52Ia7UbKwJAPWaZNA7\nNUH32VxoTkQRhcmMiJTAZEZESmAyIyIl8AWAQdoJsbvn3y310Rb83y09I/W5f0C6EI97aa/UZ9dj\nwTnW3gztj0NUVOCTVgG5YBxtYvIrEcAXAEQUYZjMiEgJTGZEpATWzAxKGrdEiF27Fkp9tDUybX3M\nqKoL4kLzUE6G1KMdCxC68bz9+UmpbeLg7iF5FrUerJkRUURhMiMiJTCZEZESmMyISAnKvQAIVfG8\n+LNvhNjfDrEUeYy8GHJpdh5JCtKuI/f9p3iQ0HuPjrhBz8Acr74otXXr1E6Iw7mbL18AEFFEYTIj\nIiUwmRGREpSrmYWLvx1iT1bVC7HZBePamkXPlHipz7l6sQ4T7Wcdt3ZH2KNn6qQ+GemJUhtZk5Gf\nC60zLrcQpye1u0HPwPzRzyYK/6aZJP7FN7VSn7sm/bsQu5xFus9izYyIIgqTGREpgcmMiJTAZEZE\nSuALAINq6sTJuJ0T9SfjGtlxwt8ERSOF3Tq3eKzd0Yp6qc+Cdw4J8a458snz2h1gm7xXhPjDr2qk\na+7pZ9cd35enLgjxxNXyjrqHl+UJcXml/GcwchQZhV6wdg7+pqZBiHt0jtO9hi8AiCiiMJkRkRKY\nzIhICZapmX11VqyX9EmTayVPv39UiH99b0azxhVpJq7/VGp7++HhN71GO+kSAAY8+qYQf7vpYd1n\nH6uQJ+x26RgrxInt2kp9rO7vX7uEeGivpBYaiXnaOmk4F5EbwZoZEUUUJjMiUgKTGREpgcmMiJRg\nmRcAkUq7+wUA2DRbYLTGwnikyJz3nhAfWXFfC43Ev6SsJ4TY9eFS3Wsam65IbW1jxM89I/7jL1If\n55P36N5bu3OLkV1b+AKAiCIKkxkRKcHQB7jly5fjwIEDaGpqwiOPPIKBAwdi/vz58Hq9SElJwYoV\nK2CzBeeQBiIiM3RrZvv27cO6deuwdu1auFwuPPDAAxgxYgSysrKQm5uLoqIipKWlweFw3PAeZmpm\n5y/KtSTtzqkq8PfnLKsUF58P+465iZgFr+0X4k0PDdO95vb5O4R4x7y7pD5GFsJ7Gr1CHNs2Wvca\n7WJ+wNiCfjM7sNKNHTl9QWrTrDNH/27twzSaINbMhg8fjlWrVgEA2rdvD7fbjZKSEuTk5AAAsrOz\n4XQ6b3YLIqKQ001m0dHRiIu7uk3H5s2bkZWVBbfbfe1rZXJyMqqrq0M7SiIiHYZfAHzwwQfYvHkz\nnn76aaE9DDM7iIh0GUpmH330EdasWYO1a9ciMTERcXFx8Hg8AIDKykrY7fqb9RERhZLuC4C6ujo4\nHA689tprSE5OBgA89dRTGDZsGMaPH4/FixcjIyMDEydOvOE96i9Ze1W+1eS88KEQ7/5Flu412mI/\nIBf8Pz9+XuozuGdHIX7n0Gkh3nXsW+ma3/xwoO54kr43R4gTB4+S+nzzyo9070Nk9AWAbrcdO3bA\n5XLh8ccfv9b23HPPYdGiRSguLkZ6ejry8/NND5SIKBh0k9mkSZMwadIkqX3Dhg0hGRARkRn8vkdE\nSrDMQvMHNTWfN/1M8Ey6/wUhdr37i2aNqzkaLsl/qLhbDH65v4524S0A3KbZZddfPUxvh1hArpFp\n62NGaHcABoBJv/1EiH/3kDwW7UTfr6vkU6h62Tm5lfRxoTkRRRQmMyJSApMZESmByYyIlBCWFwBb\nD1YIcV7/LqF+JIXRqyVfS23TvtcrbM838vJIq7LWI8SpHWJv0FN9738p/n7O2/i51OfQ0txwDUfC\nFwBEFFGYzIhICUxmRKQEy0yaVYH2rzIqKuoGPW+uyRv4wnztDrEAsHTKHUI8fmDXgMey+INjUtui\nsX2FOE8ziRYAKirEycCfPfuDgJ9t1KytpUL88oQBIXsWhR9rZkQUUZjMiEgJTGZEpATL1MxqGxqF\nuEOctU/xNnJ6VNUF+bQhe3v904b+8t9VQtwvVT71OT2pnRBrTygCgBUf/kOIHx8lz/3qo1nUrq2R\naetjRmnrWG+8+r7Ux7Vjrql7U2RhzYyIIgqTGREpgcmMiJTAZEZESmiRFwA1dXJhvE/uU0Ls2rtc\n6uO+7BXidrZoqc9nmt1Vh5jYXdVqOhWsl9q+3fRwwPcZ+vQuqe3vvx5nakyB0i7sBoBGr/ij161T\nO6lPuPgbXyQvPrcSvgAgoojCZEZESmAyIyIlWGbSLIXe/n+4pLbPq8Ua49YD4kZ9O2aODOmYrrev\n/JzUlpEmThie8dZBqU87m1hUebXgDqnP/uPinz23cKsQu96ebnicFF6smRFRRGEyIyIlMJkRkRKY\nzIhICXwBEETeK+JfZXQbczvNfnnqghDf3q297jWeRq/U1mX0L4XYVbJK9z5DntopxKOGyLvTmtnJ\n1d/Lh2HfSQr4PlZ30c8Pe7zRCjb5xRcARBRRmMyISAlMZkSkBH6ZDyKzNTKtiav3CvEfHh8j9enb\nRZxMGttWXnSfOHiUEH9dJe9G28seL8TaU5SS8lZK1zyVc5sQ+1uQra2RqVgf84f1sear0iz675Fs\nbME/P5kRkRKYzIhICUxmRKQE3S/4brcbhYWFOHfuHC5duoSZM2eiX79+mD9/PrxeL1JSUrBixQrY\nbDa9WxERhYzupNkdO3bg9OnTmD59Ok6fPo2HH34YQ4cORVZWFnJzc1FUVIS0tDQ4HI4b3sPIpNmz\n58WiX1pHc7t8rt5bLsSzR/c2dR8rqXM3Sm2J7VruKL5T37qFOFg7xD72h1KprXdn8d5zxpj773ms\nok6ItS9Qcn/ziXTNn34enB1DLjddEWJbDL8QBSJok2bz8vIwffrV7VEqKiqQmpqKkpIS5OTkAACy\ns7PhdDrNj5SIKAgMv0cuKCjA2bNnsWbNGvzkJz+59rUyOTkZ1dXVIRsgEZERhpPZpk2bcOTIEcyb\nNw/XfzMNw9JOIiJdusmstLQUycnJ6NKlCzIzM+H1ehEfHw+Px4PY2FhUVlbCbrc3eyBma2RaVqqR\nvf35Salt4uDuuteVV9YLce/UBN1r/J141TnxFt3rzNDWyPztEHtn72Td+2hrZC89IC9g195buwgf\nMLYQX1sj0wpWfcyfYNXItLXTlqybWpHu3/L+/fuxfv3Vo85qamrQ0NCAkSNHYufOq7sr7Nq1C2PG\nyDPUiYjCSfeTWUFBAZ588kk4HA54PB48/fTTGDBgABYsWIDi4mKkp6cjPz8/HGMlIroh3WQWGxuL\n559/XmrfsGFDSAZERGQGJ7wQkRJaZKdZ7SRCILwTCbWr8u1+dn24otk1tk2QdsRQkeviZaktKV5/\nRciqj8QJzt9L7yj10b5IOHJafgGQ2VX/BcA9RR8K8dwf9BHivP5ddO9BLYM7zRJRRGEyIyIlMJkR\nkRJ4OhM126QNn0ptxT8ZHvB9/E2IjdKUKv3Vx26d8bYQZ9wu17+S24t10TcfGhbw+MJJW7MFrF+3\nPVcv1k6TE4Kzkw5rZkQUUZjMiEgJTGZEpAQmMyJSAs/FUsTxavkYuZ4p8X56ih58bb8QawvlL0+Q\nd7LQamcLzo+Rkd0vtMV+ADixZqIQ/3zLIanPohxxkmzhe0eE+Ln7MqVrPI1eIfZ3nF+oWK3YX9sg\n7tjRIU7esSNYBX+z+MmMiJTAZEZESmAyIyIlcNLsDTR5xcXwMdHM+zei/bsCQvf3Ne6lvVJbn64d\nhPg3Pxwo9dn/D5cQD/tOUnAHppCGS/Iv7AW32BasnaGN4KRZIoooTGZEpAQmMyJSApMZESkhLC8A\njtdodnZtLx+BZmaSoL+hR2m2WfBXnNZicb95zB41d6yiToj9HQen3SE2tVOc1KcoX5zYe/pbt9RH\nr+C/cf8JqW3KsFtveg2FB18AEFFEYTIjIiUwmRGREjhpNsJVak6qSvVzUpWepIlrpTbX29NNj+l6\nOw5XCLG/U5SMLBqnq07UNEhtt3aW65BWwpoZEUUUJjMiUgKTGREpoUVqZn89WiX1yc6wC7GRzQbr\n/RTjEox+wSbTglFnA4Dc33wixH/6+UjTY9KjnUdmZg6Zdl4c4H9uXEtZvbdcaps9urcQv/ih3OcB\nzWlWRmpoF/387lXXXRJiI5uDGsGaGRFFFCYzIlICkxkRKYHJjIiUEJYXAO5G8RHaxeBGvXPotBCP\nH9jV9JhIHS15ilJL0v65L2hOUAKAFz4+LsRL8/qFckghwRcARBRRmMyISAmGkpnH48HYsWOxdetW\nVFRUYMqUKXA4HJgzZw4uX74c6jESEekyVDN74YUXsHfvXvz4xz/Gp59+iqysLOTm5qKoqAhpaWlw\nOBw3vV47v27QwvelPl8suTewkQP47Ph5qW1Iz45C3Ngkb87YNoYfSEkUrAmxZ1zyxpDpSe2E2PH/\nDkh9/uv/fDfgZ2kV/f8yqe2Xd93W7Pu2tKDVzMrLy1FWVoa7774bAFBSUoKcnBwAQHZ2NpxOp+lB\nEhEFi24yW7ZsGQoLC6/FbrcbNpsNAJCcnIzq6urQjY6IyKCbJrNt27Zh8ODB6N69u99/H4ZZHURE\nhtz02+iePXtw8uRJ7NmzB2fPnoXNZkNcXBw8Hg9iY2NRWVkJu91+s1sQEYWF4Umzq1evRteuXfHZ\nZ59h2LBhGD9+PBYvXoyMjAxMnDjxpteGc6fZx7cdFuIX8/uH7+EGJI1bIsSuXQt1r/n71y6p7ce/\nFXecOLLivuYNjELm9QPijh2Tvyvv2HFZ86LKZuAl1dYvTgnxhEHdpD47vzwrxD+4PU33vlYTskmz\ns2fPxrZt2+BwOHD+/Hnk5+cHegsioqAzvPnX7Nmzr/3zhg0bQjIYIiKzOOGKiJRg2dOZ7ly8W2rb\nPnuUEA+Y9ZbUp2rj1MAfRkGn3Yk0XoEdgI1MiLW6qguXpDZ7+1taYCTGcaE5EUUUJjMiUgKTGREp\ngcmMiJRg2RcAZu39qkaIR/fpHL6HG+C6KG6ZlBRva6GRtDwzE0WtxsiEWCOMvDDxXBZ3lv2vz08K\n8cP/q6epZ5sRzpchfAFARBGFyYyIlMBkRkRKaNU1s4/LaqS2UbdZq0ZG6gjVDrHBkvS/X5LaXNsf\nC8mzKs57pLYuHWND8izWzIgoojCZEZESmMyISAlMZkSkBMu+ADhefVFq65kSH4TRAEn3PS/Ervf+\nb1DuS5HHyMRfMzuIaHeIBYCtpVVCnJIoTrhenNtP976tEV8AEFFEYTIjIiUwmRGREixbMwuWPUfl\nQ4qTYtsK8R23dgzXcKiZrlyRf1zbtIlqgZEEl5FTlMor64W4d2pCwM85eqZOastITwz4Pv7UNjQK\ncYe4tjfoeXOvOL8W4jljehm6jp/MiEgJTGZEpAQmMyJSApMZESlB+RcAoeLvKLx9i3ICvk8oJwe3\nNnXuRqktsZ25IrKVaHeIjbVFS30eeesLIZ5/13ekPtqCf5/H3xHifYtzpWuSE/R3Mm7UTPxta2DH\nX22xH5AL/vPfPSL1WX5/pu69tThplogiCpMZESmByYyIlMCamUlmj7m/7z+dQvzeoyOCNiYKLk+j\nptbVVq51bf3ilBBPGNRN6rP+b8eF2N8pSov+9N9C7G/RuLZG9tWL46U+ZpSerBXiAd07BOW+wcKa\nGRFFFCYzIlICkxkRKaFFamaTfy+fcvP61OCccnPxkmYjvFsMfuGmiLZ6b7nUNmlgVyF+/fNTUp9f\n3nVbwM8ye4rSufrLQmxkDpkKWDMjoojCZEZESmAyIyIl6H4bLSkpwZw5c9CnTx8AQN++fTFt2jTM\nnz8fXq8XKSkpWLFiBWy2yPj+TkTWpPsCoKSkBG+88QZeeulfRcsnnngCWVlZyM3NRVFREdLS0uBw\nOG54DxUnzQbLGZdbaktPaqd7XVLWE0Ls+nBp0MZEVz2xQ5zIujRPnshqZIdYK9FOkAXkSbKHvpH7\nDOwh9vm6St4goZc9NBskhPQFQElJCXJyru4QkZ2dDafTqXMFEVFoGcp5ZWVlmDFjBmprazFr1iy4\n3e5rXyuTk5NRXS3vs09EFE66yaxnz56YNWsWcnNzcfLkSUydOhVe77/WrIVhmhoRkS7dZJaamoq8\nvDwAQI8ePdC5c2ccOnQIHo8HsbGxqKyshN1uD+ihJ2oapLZbO8cFdA+jzNakwsXsWFgjax7tz+Af\nvqyQ+virkWlpa2RmNyDQCtYpSmYWkRs5VclffSznhQ+FePcvsnTv4+/DUFSUudO2dGtm27dvx7p1\n6wAA1dXVOHfuHCZMmICdO3cCAHbt2oUxY8aYejgRUbDofjK75557MHfuXOzevRuNjY341a9+hczM\nTCxYsADFxcVIT09Hfn5+OMZKRHRDusksISEBa9askdo3bNgQkgEREZnBFQBEpATL7jR79rxHakvr\nGBuE0bQ+Td4rUltMdHD+P/S+pvB97+1dgnLf1iacL6XM0p6idLRCfEnQ0jvEal+2GXm5VVUr/57b\nO4i/59w1g4giCpMZESmByYyIlGDZmlmwfOOnFtKjBWsh2hPMzZ5erq2flFXWS30yu7bXvc/AJ/4k\nxIeWyqditzYNl+QfuLgI3XFYu2jc34RYI78PX2juM6hH+OpzrJkRUURhMiMiJTCZEZESmMyISAnK\nvwAIFe8V+a8tuo3+an/tBFh/k1+19zZy38OnLkht/bvpvwBQQW1DoxC7L3ulPpE64doMbbEfCF7B\n/+WP/yHEF/wkh4U5fYWYLwCIKKIwmRGREpjMiEgJrJmZ9Mt3vpTaisbfHvB9/lh6Rmr7twHputeN\n+I+/CLHzyXsCfjb9y0U/P6TxRos1OrQLsP3tpNrFQE1PWxs0siOs9hSlUJ2g5E/Sfc9LbX3vHCzE\na6Z8V+ozpGdHIWbNjIgiCpMZESmByYyIlMBkRkRKaJEXAAdPnJf63HFrR6ktGL6/aq/U9uc5o0Py\nLAq9c/WXpbbkBFuz76vdzQQwv6NJMGiL/YCxgr8e7XFwALDxoeFCHKyjGMv97OTSOzVBiB/7Q6nU\np7vmZchT37/N0PP4yYyIlMBkRkRKYDIjIiVw0qxJlX5OlUntEPhi5lAu6iW1zH/3iBAvvz+zhUbi\nf0Ls/vWPCLG2PgYAzrJzQhzXNlrqo90VeGA3+T7+8JMZESmByYyIlMBkRkRKYDIjIiVE5vlbQVDv\n561Gqom6/V2T/l1qO/jHxULs7yiwo2fqhDgjPTHwh1Or0pIFf+0OsdrdLwC54K8t9gPAiNuShbji\nvPwizcgOIv7wkxkRKYHJjIiUwGRGRErgpFmiCOLv1736wiUhtvuZ/L1k9zEhzu1tl/ps+OyUEP90\naDepj/be/upjhe+Jk4NfHG+sVshPZkSkBCYzIlICkxkRKcHQPLPt27fj1VdfRUxMDB577DFkZGRg\n/vz58Hq9SElJwYoVK2CzNX+DPCIis3RfALhcLhQUFGDLli1oaGjA6tWr0dTUhKysLOTm5qKoqAhp\naWlwOBw3vAdfAFCVn11G/BWa6apXnF9LbY+M6NUCIzFuxV/LhDg/M03q0ydNnFirLfYDwHP3iQX/\noB0153Q6MWLECCQkJMBut+PZZ59FSUkJcnJyAADZ2dlwOp3GnkZEFCK6Oe/UqVPweDyYMWMGLly4\ngNmzZ8Ptdl/7WpmcnIzq6uqQD5SI6GYMfYA7f/48Xn75ZZw5cwZTp04V5qqEYZoaEZEu3WSWnJyM\nIUOGICYmBj169EB8fDyio6Ph8XgQGxuLyspK2O3yBLrrBemUe2rFeiSzPhaIOWOsXR/zx+gpStcz\nOiHWCN2a2ejRo7Fv3z5cuXIFLpcLDQ0NGDlyJHbu3AkA2LVrF8aMGRO0ARERmWFoOdOmTZuwefNm\nAMCjjz6KgQMHYsGCBbh06RLS09OxdOlStG3b/DP9iIjMCsvaTCKiUOMKACJSApMZESmByYyIlMBk\nRkRKCNsMsCVLluDgwYOIiorCwoULMWjQoHA9OiDHjh3DzJkz8dBDD2Hy5MmoqKiw/KL65cuX48CB\nA2hqasIjjzyCgQMHWnbMbrcbhYWFOHfuHC5duoSZM2eiX79+lh3v9TweD+6//37MnDkTI0aMsOyY\nS0pKMGfOHPTp0wcA0LdvX0ybNs2y4/2nZm9o4QuDkpIS389+9jOfz+fzlZWV+X70ox+F47EBu3jx\nom/y5Mm+RYsW+TZu3Ojz+Xy+wsJC344dO3w+n8/3/PPP+954442WHKLE6XT6pk2b5vP5fL5vv/3W\nd9ddd1l6zO+9957vd7/7nc/n8/lOnTrlGzdunKXHe72ioiLfhAkTfFu2bLH0mPft2+ebPXu20Gbl\n8fp8V392x40b56urq/NVVlb6Fi1aFPCYw/I10+l0YuzYsQCA3r17o7a2FvX19eF4dEBsNhvWrl0r\nrGiw+qL64cOHY9WqVQCA9u3bw+12W3rMeXl5mD59OgCgoqICqamplh7vP5WXl6OsrAx33303AOv/\nXGhZfbzB2NAiLMmspqYGSUlJ1+JOnTpZcnF6TEwMYmPFZTdWX1QfHR2NuLir52pu3rwZWVlZlh8z\nABQUFGDu3LlYuHBhqxjvsmXLUFhYeC22+pjLysowY8YMPPjgg/j4448tP97rN7RwOBxwOp0Bj7lF\nVk36Wuk8XSuP+4MPPsDmzZuxfv16jBs37lq7Vce8adMmHDlyBPPmzbP8xgXbtm3D4MGD0b17d7//\n3mpj7tmzJ2bNmoXc3FycPHkcCC91AAABtElEQVQSU6dOhdfrvfbvrTbef2ruhhZhSWZ2ux01NTXX\n4qqqKqSkpITj0c0WFxcX0KL6lvDRRx9hzZo1ePXVV5GYmGjpMZeWliI5ORldunRBZmYmvF4v4uPj\nLTteANizZw9OnjyJPXv24OzZs7DZbJb+O05NTUVeXh4AoEePHujcuTMOHTpk2fECwdnQIixfM0eN\nGnVtYfrhw4dht9uRkJCgc5U1WH1RfV1dHZYvX45XXnkFHTt2BGDtMe/fvx/r168HcLX80Bo2Lnjx\nxRexZcsWvPXWW5g4cSJmzpxp6TFv374d69atAwBUV1fj3LlzmDBhgmXHCwRnQ4uwrc1cuXIl9u/f\nj6ioKDzzzDPo169fOB4bkNLSUixbtgynT59GTEwMUlNTsXLlShQWFlp2UX1xcTFWr16NXr3+tWXM\nc889h0WLFllyzB6PB08++SQqKirg8Xgwa9YsDBgwoNVsXLB69Wp07doVo0ePtuyY6+vrMXfuXFy4\ncAGNjY2YNWsWMjMzLTvef2ruhhZcaE5ESuAKACJSApMZESmByYyIlMBkRkRKYDIjIiUwmRGREpjM\niEgJ/wMIUJcAnJZPVgAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 360x360 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "metadata": {
        "id": "qWaTHLTH5rGw",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## Look at top co-evolving residue pairs"
      ]
    },
    {
      "metadata": {
        "id": "9A2yeOJ8uPNM",
        "colab_type": "code",
        "outputId": "98bd9e49-488e-40e2-f575-cd4ee7a32395",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 359
        }
      },
      "cell_type": "code",
      "source": [
        "######################################################################################\n",
        "# WARNING - WARNING - WARNING\n",
        "######################################################################################\n",
        "# - the i,j index starts at 0 (zero)\n",
        "# - the \"first\" position = 0\n",
        "# - often in biology first position of a sequence is 1\n",
        "#   for this index use i_aa and j_aa!\n",
        "\n",
        "# adding amino acid to index\n",
        "mtx[\"i_aa\"] = np.array([alphabet[msa['msa_ori'][0][i]]+\"_\"+str(i+1) for i in mtx[\"i\"]])\n",
        "mtx[\"j_aa\"] = np.array([alphabet[msa['msa_ori'][0][j]]+\"_\"+str(j+1) for j in mtx[\"j\"]])\n",
        "\n",
        "# load mtx into pandas dataframe\n",
        "pd_mtx = pd.DataFrame(mtx,columns=[\"i\",\"j\",\"apc\",\"zscore\",\"i_aa\",\"j_aa\"])\n",
        "\n",
        "# get contacts with sequence seperation > 5\n",
        "# sort by zscore, show top 10\n",
        "top = pd_mtx.loc[pd_mtx['j'] - pd_mtx['i'] > 5].sort_values(\"zscore\",ascending=False)\n",
        "top.head(10)"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>i</th>\n",
              "      <th>j</th>\n",
              "      <th>apc</th>\n",
              "      <th>zscore</th>\n",
              "      <th>i_aa</th>\n",
              "      <th>j_aa</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>1056</th>\n",
              "      <td>21</td>\n",
              "      <td>28</td>\n",
              "      <td>0.876234</td>\n",
              "      <td>2.977138</td>\n",
              "      <td>K_22</td>\n",
              "      <td>E_29</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>881</th>\n",
              "      <td>16</td>\n",
              "      <td>58</td>\n",
              "      <td>0.846978</td>\n",
              "      <td>2.960960</td>\n",
              "      <td>R_17</td>\n",
              "      <td>A_59</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>142</th>\n",
              "      <td>2</td>\n",
              "      <td>26</td>\n",
              "      <td>0.817712</td>\n",
              "      <td>2.943404</td>\n",
              "      <td>A_3</td>\n",
              "      <td>M_27</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>999</th>\n",
              "      <td>19</td>\n",
              "      <td>50</td>\n",
              "      <td>0.784957</td>\n",
              "      <td>2.921954</td>\n",
              "      <td>I_20</td>\n",
              "      <td>G_51</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>820</th>\n",
              "      <td>15</td>\n",
              "      <td>41</td>\n",
              "      <td>0.704959</td>\n",
              "      <td>2.860163</td>\n",
              "      <td>K_16</td>\n",
              "      <td>I_42</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>488</th>\n",
              "      <td>8</td>\n",
              "      <td>45</td>\n",
              "      <td>0.663213</td>\n",
              "      <td>2.821590</td>\n",
              "      <td>E_9</td>\n",
              "      <td>P_46</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>251</th>\n",
              "      <td>4</td>\n",
              "      <td>22</td>\n",
              "      <td>0.658282</td>\n",
              "      <td>2.816702</td>\n",
              "      <td>I_5</td>\n",
              "      <td>V_23</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1278</th>\n",
              "      <td>27</td>\n",
              "      <td>37</td>\n",
              "      <td>0.599291</td>\n",
              "      <td>2.752054</td>\n",
              "      <td>M_28</td>\n",
              "      <td>V_38</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>934</th>\n",
              "      <td>18</td>\n",
              "      <td>26</td>\n",
              "      <td>0.583831</td>\n",
              "      <td>2.733035</td>\n",
              "      <td>V_19</td>\n",
              "      <td>M_27</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>326</th>\n",
              "      <td>5</td>\n",
              "      <td>42</td>\n",
              "      <td>0.553874</td>\n",
              "      <td>2.693391</td>\n",
              "      <td>Y_6</td>\n",
              "      <td>H_43</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "       i   j       apc    zscore  i_aa  j_aa\n",
              "1056  21  28  0.876234  2.977138  K_22  E_29\n",
              "881   16  58  0.846978  2.960960  R_17  A_59\n",
              "142    2  26  0.817712  2.943404   A_3  M_27\n",
              "999   19  50  0.784957  2.921954  I_20  G_51\n",
              "820   15  41  0.704959  2.860163  K_16  I_42\n",
              "488    8  45  0.663213  2.821590   E_9  P_46\n",
              "251    4  22  0.658282  2.816702   I_5  V_23\n",
              "1278  27  37  0.599291  2.752054  M_28  V_38\n",
              "934   18  26  0.583831  2.733035  V_19  M_27\n",
              "326    5  42  0.553874  2.693391   Y_6  H_43"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 12
        }
      ]
    },
    {
      "metadata": {
        "id": "wqG93dC12CKx",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## Explore the MRF"
      ]
    },
    {
      "metadata": {
        "id": "k6BsheyNx3ID",
        "colab_type": "code",
        "outputId": "79122395-0c55-49aa-a00e-2eb2331a1fc2",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 333
        }
      },
      "cell_type": "code",
      "source": [
        "def plot_v(mrf):  \n",
        "  al_a = list(alphabet)\n",
        "  v = mrf[\"v\"].T\n",
        "  mx = np.max((v.max(),np.abs(v.min())))\n",
        "  plt.figure(figsize=(v.shape[1]/4,states/4))\n",
        "  plt.imshow(-v,cmap='bwr',vmin=-mx,vmax=mx)\n",
        "  plt.xticks(np.arange(v.shape[1]))\n",
        "  plt.yticks(np.arange(0,states))\n",
        "  plt.grid(False)\n",
        "  ax = plt.gca()\n",
        "  ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x,y: mrf[\"v_idx\"][x])) \n",
        "  ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x,y: al_a[x]))\n",
        "  \n",
        "plot_v(mrf)"
      ],
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1sAAAE8CAYAAAA2UqvOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3WlgVeXZ9v0jJEwhzBAGEwhQJmUS\nBASZRBlE0ApaLYJFxbs4gC+9QZBSRagKOFRBizjVCrWogIoMokhBBIyCTFFBGUNE5jEMAsl6PvjA\ng+nae6/rlEXfpP/fpybl8LjYOffa+2KtrB3neZ4nAAAAAMB5Veg/vQAAAAAAKIjYbAEAAABACNhs\nAQAAAEAI2GwBAAAAQAjYbAEAAABACBJCb8jIsOXKlnWOzPjsIlNVr44HnDOzl7qvT5LatjXFtH27\ne2bdOltXdrZ7pmvXC9d1+LCty2rZMvdMSoqtq5Dhnz+sj8fJk7ZczZrumblzbV1pae6ZgwdtXe3a\n2XJffumeqVvX1vX+++4Z6ywecD8sSpJuv909s22brev7790zVavauvbvd88kJ9u6Nm2y5YoWdc8k\nJdm6fltnpXuoTh1T14YdJU25UqXcM2vXmqrUpfUR58y6rba/l/XYvXGjeyY319ZlmaseacY3LevX\nO0dun3uTqWrECFNML7/snune3dZVubJ7plgxW1dcnC23fLl7pn59W5fl9Xb1altX8+b+3+fMFgAA\nAACEgM0WAAAAAISAzRYAAAAAhMB5szV79mxdcskl2m+5gB0AAAAA/kuYNlupqamaP39+GOsBAAAA\ngALBabN18OBBrV27VsOHD9ecOXPCWhMAAAAA5HtOm60PPvhAHTp0UNu2bbV161bt2rUrrHUBAAAA\nQL7mtNmaPXu2unfvrvj4eHXt2lVzrR+eAwAAAAAFXOAPNd65c6fWrFmjsWPHKi4uTidOnFDJkiV1\nu+UTLAEAAACggAu82Zo9e7ZuvfVWDR8+XJLkeZ46d+6szMxMVatWLbQFAgAAAEB+FPgywjlz5qhn\nz55nv46Li9Ovf/1rbpQBAAAAAD4Cn9l65513/u17995773ldDAAAAAAUFM6fswUAAAAAiC3O8zwv\n1IaMDFuufPkL1jXjcCfnTK/kJaau7yq3NeVqL5zsnDl0y+9NXZmZ7pmGm98zdQ1bdr1zZvRoU5Wy\nsmy506fdM8eO2boWLHDPFCtm66pRw5a74gr3TNlZfzd1TT7xO+fMRReZqtS9tO05rdxc58ibO9ub\nqm6+7rh7aNo0U5f69DHFXnqtsHPmrtt+NHXpxRedI7t+M9BUVWnas86ZxU3uN3UlJppiSktzz1g/\nwaXBjyudM5M+b2bqurtpuimnhQvdMy1amKq+qXqVc6Z+5nxT16rkLqac5bVs7VpTle7s7X6smv9J\ncVPX5Ze7Z0ovmGHqWvOrXqZc443ufVktbV2bNrlnWrc2VWn/fluuUrlTzpkjJ9xfWyRp8WL3TMOG\npipVr+7/fc5sAQAAAEAI2GwBAAAAQAjYbAEAAABACALfjTArK0s9evRQgwYNJEknT55UnTp1NGrU\nKMXHx4e2QAAAAADIj5zObNWoUUNTpkzRlClT9Oabb+rUqVN6//33w1obAAAAAORbv+gywkaNGmnb\ntm3nay0AAAAAUGCYN1unTp3Sxx9/rEsuueR8rgcAAAAACoTAv7MlSVu2bFHfvn0lSRs2bFD//v11\n9dVXh7IwAAAAAMjPnDZbZ35nS5IGDRqkGtZPRQUAAACAAs58GeHQoUP15JNP6vhx908JBwAAAICC\nzrzZSk1NVZcuXTRp0qTzuR4AAAAAKBACX0aYkpKimTNn/ux7f/jDH877ggAAAACgIHD6na0Lac3u\nKs6ZxsowdVWqZAh5piqtXWvL1e7XzzlTWj+ausqWLeoeuvhiU9eQ1u6ZrCxTlWql2B6PrD3uj0fd\nyodMXZf2znYPVa5s6tq01fZh5GUXzjDlLO64wz2zYoWx7LQt9ml8e+dMnTq2Lu3fbwwanLY9IE2b\nFnbOfLHWcMyR1NyQ+fBDU5V+fcf9zpkym21djd99xJTb1Odh50y87TBgcnc340fF7DQWDhrknvns\nM1NVZqZ7Zl+JLqauNl++Ysrtue5O50zNmqYqrfm2uHPG8hhKUpfDbztnsq+5ydRVz/3wJkn659e9\nnDPFvrB13dBki3Mmp5DtHgyVTm435ZSQ4hxJTLRVXXute2b3bltXJL/oc7YAAAAAAP7YbAEAAABA\nCNhsAQAAAEAIAm+2srKyVL9+fa1fv/7s92bOnPlvN80AAAAAADie2frVr36lp556Kqy1AAAAAECB\n4bTZuuSSS5SYmKjly5eHtR4AAAAAKBCcf2dr8ODBeuaZZ+R5xnufAwAAAMB/AefNVlpami6++GLN\nnTs3jPUAAAAAQIFguhvhvffeqxdffFGnjR96CQAAAAAFnWmzVaFCBV199dWaNm3a+V4PAAAAABQI\n5s/ZuuOOO7Rz587zuRYAAAAAKDASgv7BlJQUjR079uzXJUqU0LJly0JZFAAAAADkd+YzWwAAAACA\nyAKf2brQGn8wzj3UtKmpq3hxQ+iYqUopKbacXnvNPZOVZao68JsxzpkjCbVNXb8q456ZPNlUpVat\nippyO3a4Z7p2LW3qqrhqiXNmU/3upq5jxhm+kDIz3TOff27ramU7fKhNzmLnzDtb25u6Lk02xWwm\nTDDFdjca5py5/HJTlVSvn3Mke6qtqmTuIefM11/bjgOVBzxsyr3yrHvm5ptNVSbbC1U35VJl/JUF\nywy3aGGq+tWv3DOWl3VJapNmy336qXvG+LZKH37onjG9FzNKmj/DlNvXoZcpV6+ee2bzZlOVVMj9\nPEr85u9sXVu3mmJr9qc6Z9auNVWpd2/3zLff2roqVfL/Pme2AAAAACAEbLYAAAAAIARstgAAAAAg\nBIF/ZysrK0s9evRQgwYN5Hme4uPjNWDAALVq1SrM9QEAAABAvuR0g4waNWpoypQpkqTMzEwNGDBA\nTz/9tOpZfvMPAAAAAAow82WE1apV04ABA/TGG2+cz/UAAAAAQIHwi35nq0GDBtq4ceP5WgsAAAAA\nFBi/aLN19OhRxcfHn6+1AAAAAECB8Ys2WxkZGapfv/75WgsAAAAAFBjmzVZmZqZee+019evX7zwu\nBwAAAAAKBqe7EW7ZskV9+/bVyZMnlZOTo4ceekhVq1YNa20AAAAAkG8F3mylpKRo1apVYa4FAAAA\nAAqMOM/zvFAbMjJsuTlz3DNNm9q6srPdMxUq2LrmzbPlHnnEPTN6tK0rOdk9U62arSsz0z1zzz22\nLmuubVv3zJEjtq7UVPdMWpqtq3BhU+z9b+s6Z6wfxVf7zT+7h8qUsZU1bmyKbUlxn48ar/7J1PXD\nPWOcM1lZpiodO2bLtV/2uHuoZ09TV/pB91ls2eCoqeuVaSWcM3fuf8LUZX7wLcduyzFYkm680T0z\na5atq1s3U+xU05bOmcKffGzqUmKie+bwYVPVocu7mHKlT+xyD734oqlL/fu7Z6ZONVUdv2+ocyYn\nx1SlpK/SbUHLgTgpydb16afumTvusHWVK2eKHcgt7ZwpO8HwPliSduxwz3TubOvq1cv327/oBhkA\nAAAAAH9stgAAAAAgBGy2AAAAACAETncj3Lp1qx577DHt379fubm5uvTSSzVs2DAVKVIkrPUBAAAA\nQL4U+MxWTk6OBg4cqP79+2v69OmaMWOGJOn5558PbXEAAAAAkF8FPrO1dOlS1axZUy1atJAkxcXF\naejQoSpUiCsRAQAAACCvwJutzZs3q379+j/7XrFixc77ggAAAACgIAh8WiouLk451g8mAAAAAID/\nMoE3WzVr1tS6det+9r2TJ0/q22+/Pe+LAgAAAID8LvBm64orrtD333+vhQsXSpJyc3P1xBNPaO7c\nuaEtDgAAAADyq8CbrUKFCumVV17RW2+9pZ49e6p3794qWbKkBg0aFOb6AAAAACBfcvqcreTkZL3w\nwgthrQUAAAAACgzu2w4AAAAAIYjzPM8LtSEjw5arWNE9s3atqerN/Z2cMzdXXWLqWpzb1pTLzXXP\nXFnR9tj/a08D50yTJqYqk7I6YMptOVjWlPvqK/fM+++bqvTnP7tnli2zdVWtass1z5zhHvr1r01d\nM96Nd85Y/16tTtue05Yn5xeJ7U1VzZuccg9NnWrqOtTzdlPuyy/dM1e2Mfy9JMlwpUV6i4Gmqpb1\nDjlnfjhW2tRlvfFverp75v9+dKaz1N0rnTPLTzYzdbUqZPiLSTpQp6VzpuyXH5u6Pi16lXPmlHHs\nr9z8iimX3uBO50xamqlKBw+6Z44ds3VduvFt95DxM2KPd+tlyu3e7Z6pvn6+qWtebhfnTOvWpiqV\n3m57n5lT3/19pvV+fJYZPnnS1lU6wiGfM1sAAAAAEAI2WwAAAAAQAjZbAAAAABCCwHcj3L59ux59\n9FHt2bNHnufpsssu0//+7/+qaNGiYa4PAAAAAPKlQGe2cnNzNXDgQP3ud7/TjBkzNHPmTFWuXFkP\nPfRQ2OsDAAAAgHwp0GZr6dKlql69ulq1anX2e7fffrtWr16t/fv3h7Y4AAAAAMivAm22Nm/erIsv\nvvhn34uLi1Pt2rW1ZcuWUBYGAAAAAPlZoM2W53nK8fngD8/zFPbHdAEAAABAfhRos1WjRg1l5Plw\nYs/ztHHjRtWsWTOUhQEAAABAfhZos9WmTRtt2rRJixcvPvu91157TZdeeqnKlSsX2uIAAAAAIL8K\ndOv3+Ph4vfzyyxo2bJieeuopeZ6nSy+9VI888kjY6wMAAACAfCnwhxqnpqbqjTfe0KhRo1S8eHGN\nGjWKz9gCAAAAgAgCb7bOaNq0qRo1aqSePXtq3rx5YawJAAAAAPK9QJcR5jVy5MjzvY5/s0uVnDPu\niZ8kJxuDBgcP2nLXdzvlHtpqO/P40UfumapVTVWqm3zAOfPpV2VNXYULm2KaMsU907evrSs31z2z\nc6eta+tWW655iiE0daqpK1u/c8689JKpSq1ut+UsPvvMljt92n2IW8X+I77Wr7fl1q51z9SubXty\nWkbROveXXVbaObM309ZlfU5//bV7pnNnW5d2u0cqVDB2GT/Os2yS4XXTqHhx98y999q61gyy5fLc\n5yyQYsVsXdOmuWe6dbN1XUjFD+8y5TIz3d+hbi7SxdR10HD8KL37O1OXdtkejx2lGzhnZs0yVal3\nb/eM9SOEGzf2/77zmS0AAAAAQGxstgAAAAAgBGy2AAAAACAETr+zlZWVpR49eqhBg59fazlx4kSV\nKVPmvC4MAAAAAPIz5xtk1KhRQ1MsdwwAAAAAgP8iXEYIAAAAACFgswUAAAAAIXC+jHDLli3qe86H\nCNWoUUOjR48+r4sCAAAAgPyO39kCAAAAgBBwGSEAAAAAhOAXX0YoSUOHDlWjRo3O26IAAAAAIL9z\n2mylpKRo1apVYa0FAAAAAAqMOM/zvFAbMjJsuTlz3DNNm9q6li1zz3TsaOuaN8+We+QR58j2nYVN\nVdnZ7pkBA0xVWnzjROfMmnYDTV2NkzaZcptUyzlTs6apSnFzZruHunQxdVnnI/XzGc6ZxRV6mbqa\nNHHPvPuuqUq/q7nEFszNdY5kN2tvqkqKP+4emjbN1OX1u92UO3zYPVN63AhT1/IejzlnhgwxVemN\nN9wzSUm2rvLlbC/LGV/FOWfKlDFVKWXXSvdQuXK2st27bbmFC90zLVrYuma4HxezRvzVVJUy/xVT\nbmWTO50z1hm2HAes41Hry7edMwM/ucnUNfHpU6Zct+vdX2+3bzdVad2nh9xD48aZunbd734MlqQT\nJ9wzlSubqlR06wbnzIHkuqausmX9v8/vbAEAAABACNhsAQAAAEAI2GwBAAAAQAgC3yBj27Ztevzx\nx7Vv3z5JUtWqVfXwww+rnPUiWwAAAAAowAKd2crJydHAgQPVv39/vf3223r77bd1ySWX6NFHHw17\nfQAAAACQLwU6s7V06VLVrl1bl1122dnv9e/fX2HfyBAAAAAA8qtAm63Nmzerbt2f3waxUCF+3QsA\nAAAAIgm02SpUqJBOnz599uu7775b2dnZ2rlzp2bNmqXixYuHtkAAAAAAyI8CnZ6qXbu21q1bd/br\nSZMmacqUKcrJyVGu4QM9AQAAAKCgC7TZuvzyy7Vz504tPOfT2L/66isdPXpU8fHxoS0OAAAAAPKr\nQJcRxsXF6eWXX9bo0aP1/PPPq3DhwkpMTNSkSZNUrFixsNcIAAAAAPlO4M/ZKl++vJ599tkw1wIA\nAAAABQa3FAQAAACAEMR5YX9YVkaGKXakegPnTMnPPjJ1/X1HJ+fM72ouMXWtKdXWlGt88SnnzCkV\nNnVlZ7tnyn7ynq3rquudM0lv/83UdeTG2025hMDnf/+f4q9NMnUpNdU5kt2hu6kqqYTxqT9zpnPk\nSOdepqqvv3bPtPziOVOXGjc2xf620f053a2bqUqVZr/iHjJ+TMeMUrbnS8+e7pkFC0xV6rR+onNm\ncaOBpq72qw1XdvTubeqS8fL873aWdM6cOGGqUsOTK50z83Y3M3VdUy7dlDtUr6VzpvSKj01dS4pc\n5Zxpe2y+qevTEl1MuTYNDzln1mWWNnWVKuWe+fxzU5Vu0tvuIeNxcXsL22tZYqJ7pvwK23yoY0f3\nzMMPm6pOPfKYKVd463fOGe9XtU1dcbt3uYfOuQO7k4su8v02Z7YAAAAAIARstgAAAAAgBGy2AAAA\nACAEgTdbWVlZ6pnnYvyJEydq6tSp531RAAAAAJDfcWYLAAAAAELAZgsAAAAAQuB0U+stW7aob9++\nZ7/+/vvvdccdd5z3RQEAAABAfue02apRo4amTJly9uuJE90/4wQAAAAA/htwGSEAAAAAhIDNFgAA\nAACEgM0WAAAAAIQg8O9spaSkaObMmT/73sCBA8/7ggAAAACgIODMFgAAAACEIM7zPC/UhowMU+wv\nHzVwzkyYYKrSunXumSuusHWtuXaEKTe5+mPOmTp1TFW6MsP9LpN7brGd5Uxwuh/mT8q++zdT1/ar\nbzflUhdPdc5kdehj6vrnP90zgwaZqlS0iO2pP+bPcc6ZkydNVRrT4n3nzKmuPUxdTZuaYqbjzpWX\nHTF1HThd0jnz6qumKnXrZsvVL/W9c2be2otMXZbjR6d6201dk+emOmd+39rw4iJpnRqacg0/neQe\natvW1NVxkPtr9Ouvm6o0eLAtV7u2e2bVKlvXvNk57qEFC0xdd03vYsq9VO8p99CxY6Yu05PT+FFC\nby6q5Jy5OWGGqWt+Ui9Trku2oS8pydS1rqr7fKSkmKr0vfvhXpLU4BLD+4/Tp01d23cWds6kZn9j\n6lL9+r7f5swWAAAAAISAzRYAAAAAhOAXbbaysrLUs2fP87UWAAAAACgwOLMFAAAAACFgswUAAAAA\nIWCzBQAAAAAhYLMFAAAAACFgswUAAAAAIWCzBQAAAAAhYLMFAAAAACH4RZutlJQUzZw583ytBQAA\nAAAKDM5sAQAAAEAIEsIueH9LA1OuQgX3zJ13mqqUFH/cOdO/f3Fb2X2PmmJZD7ln6tQxVUlduzpH\nMjNtVUWLumfK3nabqSthtykmXXONc2TrN7aq5GT3TNHTR01dW3aUMOX+9S/3zF13map0qF0P50yx\nXFtX27a23JWFFjtnjie0N3Wt/dI9c9z98CZJKlPGllMF9yG2PvYDB7pnOrX+wNTVoYP7EB+v1tDU\ntX62KSZTW/Xqpq5773XPpLz0sKnrxhsfMeVOnnTPFC5sqpKOHXOO7GnaxVTVdLMpJt1xn3Nkyw7D\ni7SkGgnbnTMbDlYydZUr555ZV7mXqatpZVNM+sQ98kMj23wkuo+iXnvNVKUGtrf4alDf8EK92Tb4\np4vUdc5kp9Y3dSVF+D5ntgAAAAAgBGy2AAAAACAEbLYAAAAAIATOm62srCz17NkzjLUAAAAAQIHB\nmS0AAAAACAGbLQAAAAAIAZstAAAAAAgBmy0AAAAACAGbLQAAAAAIAZstAAAAAAhBgiW0ZcsW9e3b\n9+zXQ4cOVaNGjc7bogAAAAAgv3PebKWkpGjVqlVhrAUAAAAACgwuIwQAAACAEMR5nueF2pCRYcvN\nmeOeadrUVLU8qZNzptXpJaau9CJtTbmWqye7h265xdT1yvTSzpk7K7xn6voo8XrnTFqaqUrFitly\n27a5ZzZssHW9+qp75o47bF2FjP/UcnupGe6h7GxT16Rjv3PONGxoqlIbz/ac3lTV/Tm9caOpSq1b\nu2dKTv+braxJE1Psb6svdc7c/pujpi7LE2Z2jYGmqjJl3DMHDpiqVLasLVexonsmM9PW1ancSufM\nlnLNTF01dqebclq40D3TooWpKrvlVc6ZpKXzTV3KyjLFFv/qTueM9XUiOdk9Y31t//pr98ylmw2v\nY5I2NellyuXmumdqb7bNx9uHuzhnblo1wtT13e2PmXJFirhn9u83Vemo4eWlbl1bV6RjMGe2AAAA\nACAEbLYAAAAAIAROm62srCz17Nnz7NcLFizQrbfeqpMnT573hQEAAABAfma69bskbdiwQRMmTNBr\nr72mIpaLLwEAAACgADNdRrh//34NGzZMf/nLX1SuXLnzvSYAAAAAyPecN1unT5/WoEGDdM0116hW\nrVphrAkAAAAA8j3nzdaWLVt0zTXXaMaMGdq5c2cYawIAAACAfM95s1W7dm3deuut+sMf/qAhQ4Yo\nJycnjHUBAAAAQL5mvvV7165dlZqaqueff/58rgcAAAAACoRf9DlbI0eO1Ny5c5Webvx0dwAAAAAo\noJxu/Z6SkqKZM2ee/bpEiRL64IMPzvuiAAAAACC/+0VntgAAAAAA/uI8z/PCLNi2zZY7edI9s3mz\nratLZ/eH4M234kxdN1933JR79sXizpmOHU1VarhoonPmi8sHmrqa1znkHjrn7KqTNm1sOcNlspta\n9TFV1azpnol7+ilTl267zRT78+SKzpnq1U1V6rvtz+6hMmVMXZPi7zPl7r54sXvo2DFT1/EO1zhn\npk41Vem3v7Xlkg5sdw+VKmXqWrK2tHOm7fqXTF3Zv73LOZO0aY2pS5Ur23KWY2Pnzqaq9zLcP/rl\n+nYHTF3vLCprylWt6p45ccJUpfbJ37iHMjNNXfPVxZTrkmF4rejZ09SlBKcLpyRJxyukmqoKGU4b\nrF5tqlLLQl/YgpafdVKSqWpXE/f5OHjQVKW6dYxbiA8/dM906GCqOnSiqHOm9MaVpi41a+b7bc5s\nAQAAAEAI2GwBAAAAQAjYbAEAAABACAJttrKyslS3bl2tznORa69evTR8+PBQFgYAAAAA+VngM1up\nqamaPXv22a+3bdumw4cPh7IoAAAAAMjvAm+2GjdurGXLliknJ0eSNGfOHF1xxRWhLQwAAAAA8rPA\nm63ChQurcePGSv+/t8H++OOP1b59+9AWBgAAAAD5mdMHIXTt2lWzZ89WhQoVVKlSJSUmJoa1LgAA\nAADI15zuRtiqVSulp6drzpw56tLF9qF6AAAAAPDfwGmzVaRIETVv3lwzZsxQx44dw1oTAAAAAOR7\nTpcRSj9dSrh//36VLFkyjPUAAAAAQIEQaLOVkpKisWPHSpI6dOigDh06SJJatmypli1bhrY4AAAA\nAMivnC4jBAAAAAAE43wZoavqRzJMuTe/auCcKVXKVCWNHu0cubzfw6aqAyeKm3Lt2rlntm41Vamh\nIVOnjq1r++HSzplUW5W0Y4c16azWh5NswVTz387dhAmm2EU1xzhnLPMrSUr8vXvmzTdtVRfy5qot\nWphixd94xZC609SVdOh7U05Tp7pn+vc3VR07ZoqZJOUccg99+KGtzPoXS052z7z8sqmqSLvH3UMH\nD5q66tUra8plZ7tn/vpXU5XaP219A+KuRIkLViUlJdlyhrkqfscdpqpnp1VyzpQrZ6pSyw6VbcHM\nTFvOoNL+b5wziSn1TV1Z38eZcinz57uHCtnOD5W2vEFt1MjUFQlntgAAAAAgBGy2AAAAACAEbLYA\nAAAAIAROv7OVlZWlHj16qEGD//f7VPXq1dMf//jH874wAAAAAMjPnG+QUaNGDU2ZMiWMtQAAAABA\ngcFlhAAAAAAQAjZbAAAAABAC58sIt2zZor59+579unXr1rr77rvP66IAAAAAIL/jd7YAAAAAIARc\nRggAAAAAIWCzBQAAAAAhcNpspaSkaObMmWGtBQAAAAAKDM5sAQAAAEAI4jzP80ItiLPlvLHj3ENN\nm9rKli1zz3TsaOuaN8+Wq17dPZOVZar6W80xzpmKFU1V6r5lonMmvcVAU1dioimmiy92z8S/OMlW\nlprqnrn2WlvXQw/Zck2auGeys21d27e7Z8qUsXU1bmzL5ea6Z6pWtXV98ol7ppDx39R277bl+vVz\nzzz7rK2rShX3TLFitq5jx5wjv/vyflOVdTwGDXLPVJnwoK3sxhvdM7NmmarW3fiIKddw9uPuoRYt\nTF1KSnKO7KnZ0lRVcdYrptzQ9Xc6ZyxvPSTpvusNx+5p02xlaWnuGetx8de/tuX27nWO/JBbyVRV\n5dgm58zirFqmrvbVtphySklxjhzILmyq+v5794zlrZgklS7t/33ObAEAAABACNhsAQAAAEAInD5n\nKysrSz169FCDBg3keZ5Onjypu+66S506dQprfQAAAACQL/2iDzU+ePCgbrjhBrVt21bFrNfBAwAA\nAEAB9IsuIyxTpowqVqyoPXv2nK/1AAAAAECB8Is2W1lZWTp48KCqWO4IBQAAAAAFmPNlhFu2bFHf\nvn3leZ6KFi2qcePGKSHB+T8DAAAAAAXaL/qdLQAAAACAP279DgAAAAAhYLMFAAAAACFw2mylpKRo\n5syZYa0FAAAAAAoMzmwBAAAAQAjiPM/zwixYvtyWa3XZKffQokWmrk01Ozlnau1YYup6Z29bU+6G\n3ZOdM+ta/97UlZ3tnqlQwVSllBT3zOjRti5rbvVq90ypUrauut/Nds6srNLd1FWunCmmGl/OcM7s\natPL1LV1q3smN9dUpVanbc/ppz53f07/730/mromv1bUOdOokalKGzfack2auGd27rR1WY87Fpab\n7h46ZOtKS7PllhhGOCnJ1tWknMl7AAAfFElEQVSj6krnzEo1M3U1O51uyv1QraVzpsrXH5u6vql6\nlXPm5ElTlfn15fPP3TOXXGLrWrXKPZOcbOuyHOOmTrV13XabLWd5e3pzmfm2sjp1nCMr99cwVTUr\nmmHKzf++gXPmootMVSpd2j1TubKtq3Bh/+9zZgsAAAAAQsBmCwAAAABCwGYLAAAAAELgdBX6P/7x\nD7333nsqUqSITpw4oT/84Q9q3bp1WGsDAAAAgHwr8GYrKytLb731lqZPn67ChQtr69atGjlyJJst\nAAAAAPAR+DLC7Oxs/fjjjzp16qe7BKalpWmq9XYuAAAAAFDABd5s1atXT40aNdJVV12l4cOHa+7c\nuTp9+nSYawMAAACAfMvpBhnjx4/X1KlTVa9ePb388su6/fbbFfLHdAEAAABAvhR4s+V5nn788UfV\nqlVL/fr109tvv61du3Zpx44dYa4PAAAAAPKlwJut6dOn609/+tPZM1lHjhxRbm6uypcvH9riAAAA\nACC/Cnw3wp49e2rz5s266aablJiYqNOnT2vkyJEqVqxYmOsDAAAAgHwp8GYrPj5ew4YNC3MtAAAA\nAFBgON0gAwAAAAAQTOAzW1atSmaYcqsyGjhnLjU1Sd98454p0ritqevwZlNMn17ye+dMWllbV1KS\ne8Z6NWnxEwecM7m5tr/Yzp2mmJo0cc9k2MZeWU26O2cS9tm6rJ/c8EW1Xs6ZarYqVajgnqlUydb1\n3Q+25/Trr7tn+vcvauoqUsQ9s3evqUoN3A/BkqSDB205i88/d8/8vtBLpq53KtzlnLn4YlOVjh+3\n5Ro2dM9s2WLrenhWM+dMNeOBILF1S1Mux3BsrGJqkjIz3TMlSti6du+25S65xD1jfT5bjlV16ti6\nSpZ0z9SrZ+s6dsyWs/S9n9nF1FXjqHsmK8tUpWbtLjLlLKkGtWwHRq9YcedMdrapSoUL+3+fM1sA\nAAAAEAI2WwAAAAAQAqfLCMeOHauvvvpKe/bs0fHjx1WtWjWVLl1azz33XFjrAwAAAIB8yWmzNXz4\ncEnSzJkz9d1333F3QgAAAACIgMsIAQAAACAEbLYAAAAAIARstgAAAAAgBGy2AAAAACAEbLYAAAAA\nIARstgAAAAAgBE63fj+jZ8+e53sdAAAAAFCgcGYLAAAAAEIQ53meF2pDRoYtd8kl7pkFC2xdy5a5\nZzp2tHXNm2eKTb34MedMmTKmKnXfMtE9VK2araxDB/fMzJm2rpo1bbnt290zR47YulJT3TMbNti6\nDh605Zo0cc9kZ9u6LI+9dfAbN7blcnPdM8eO2bp27HDPFDL+m9ru3bac5fGwPl+qVHHPFCtm67L8\nzE6evHBdkpSc7BzZ1u1uU1X1vSudMxlFm5m6GhxNN+U2VWjpnKm19WNTlxIT3TOHD9u6srJsOcsx\n3zqLCYYLpywZSQ1fH+qcWTdqhqlLixbZcpb3OklJtq5PP3WOeKPHmKriDh8y5fTII+6Z0qVtXY0a\nuWc++MDWNXmy77c5swUAAAAAIWCzBQAAAAAhYLMFAAAAACGIutnq0aOHMjMzz37drVs3LV68+OzX\n9957r5YsWRLe6gAAAAAgn4q62WrZsqW++OILSdL+/ft1/Pjxs19L0po1a9Ssme0XXgEAAACgIIu5\n2VqxYoUk6csvv9R1112n1atXS5I2bdqklJQUJVruwAMAAAAABVzUzVbz5s21cuVPt3ddsWKFWrdu\nrZycHJ04cUJffPGFWrZ0v60qAAAAAPw3iLrZKlOmjBITE7Vr1y6tWbNGjRs3VqNGjbR69WqtWLFC\nl19++YVaJwAAAADkKzHvRtiyZUstWbJEcXFxKlasmJo1a6ZVq1Zp3bp1uvTSSy/EGgEAAAAg3wm0\n2XrzzTfVpEkTSVKzZs20aNEiVaxYUcWKFQt9gQAAAACQH8XcbDVv3lxfffXV2bsOli9fXgcPHuT3\ntQAAAAAgioRYf6BkyZL6+uuvf/a9+fPnh7YgAAAAACgIYp7ZAgAAAAC4i3lm6z/l06VxzplTCZ1M\nXVcOae0e+vJLU9eqmx4z5fp8Ptk5c6jH701d2mKLmZQqdeG6ype35bZvP7/r+P+LESNMsS8yijtn\nml/mmbqUmemeef99W9eFVLLkBav6KOV2U67T7nHneSVRPPigKfbFevfHsfnql0xdF9Qtt9hyRYo4\nR+LjbVUWx49fuC5JKmT45+SVZa4ydTVt6p6J+9B2ldCSOneacm0/f8qUszj+/7k/py0/L0m6zZD5\n58lepq7f9kkx5XTqlHvm6FFbl0HcwQOmXHbhsqZc0p/+5B6aMMHUpQoV3DN/+IOtKwLObAEAAABA\nCNhsAQAAAEAI2GwBAAAAQAgCb7b++te/6umnnz77dW5urq6//nqtX78+lIUBAAAAQH4WeLN1xx13\naP78+dq1a5ckacaMGWrcuLHq1asX2uIAAAAAIL8KvNkqVqyY7rnnHj3zzDM6fvy4Xn31Vd1///1h\nrg0AAAAA8i2n39m67rrrtGnTJo0cOVI33HCDyltvpQ0AAAAABZzTZisuLk6DBw9Wenq6+vXrF9KS\nAAAAACD/c74bYWpqqpKTk1XE8OGJAAAAAPDfglu/AwAAAEAI2GwBAAAAQAicN1spKSmaOXNmGGsB\nAAAAgAKDM1sAAAAAEII4z/O8MAvmz79wOWvXggXuma5dbV2LFtlyljUWL27r6r5lontowABT15ET\nhZ0zJaf/zdS177rbTbnyuXucM4u/rmjq2rzZPXP7bTmmLi1fboqNW9rG1mdw663umQMHbF0PPGDL\nvfCCe6b6sW9MXX96o75zJjvbVKW0NFvullvcM19+aeuqUME9U7OmrWvCBPfMIw8cNXVtyCphytWt\nesQ5s+dESVOX5TVw+nRTlbp3t+XuuMM9s3SprWvSJPdMpdW2Ny23T+tiyv3tfzPcQ6mppi799a/u\nmYQEU9WSy4c6Z6pWNVVp2TJbrm/iDOdMVstepq6PP3bPtGtnqtLp07ac5dhdtpTxvc6JE86RnGK2\nY3B8vP/3ObMFAAAAACFgswUAAAAAIQi82br55puVkfHzU9BPPfWUXn311fO+KAAAAADI7wJvtrp3\n76558+b97Hsffvihrr322vO+KAAAAADI7wJvtrp166aPPvro7NcZGRlKTk5WpUqVQlkYAAAAAORn\ngTdb5cuXV2pqqtauXStJmjdvnnr06BHawgAAAAAgP3O6QUb37t01d+5cSdLChQvVpYvt9qMAAAAA\nUNA5bbY6deqkf/3rX1q3bp3S0tJUunTpsNYFAAAAAPma02YrKSlJdevW1eTJk7mEEAAAAACicP6c\nrR49emjp0qXq2LFjGOsBAAAAgAIhwTXQqVMnrVq1Koy1AAAAAECB4XxmCwAAAAAQW5zneV6YBcOG\n2XKjR7tn1q+3dTVOO+ScWbnRdnOQZmn7TDlNn+6eycoyVW3qN8Y5UyvjPVOX2rVzjhxQWVNV2b3f\nmXJKT3fPHDli60pNdc9s2GDrOnjQlmvSxD2TnW3r2r7dPVOmjK2rcWNbLjfXPbNgga0rLc054t1x\np6kqbvw4U870ePTvb6r6YmtF50zz1S+ZunTsmHvm5MkL1yVJycnumcxMW9eNN7pnZs2ydXXrZssZ\nni+bsm2fHVrr4ErnTFalZqaulPmvmHKmY751FhOcL5yyZSTTz1mFjOcaLF2Ssuu6/6ytL9EpF7m/\nrV/+WZypq2pVU0x797pnNm60dZUq5Z45fNjWdfPN/t/nzBYAAAAAhIDNFgAAAACEgM0WAAAAAISA\nzRYAAAAAhIDNFgAAAACEgM0WAAAAAISAzRYAAAAAhIDNFgAAAACEwPgJctIbb7yhefPmqWzZspow\nYcL5XBMAAAAA5HvmzVbv3r3Vu3fv87kWAAAAACgwuIwQAAAAAELAZgsAAAAAQsBmCwAAAABCEOd5\nnvefXgQAAAAAFDSc2QIAAACAELDZAgAAAIAQsNkCAAAAgBCw2QIAAACAELDZAgAAAIAQsNkCAAAA\ngBCw2QIAAACAECT8J0ofe+wxrVmzRnFxcRoxYoQaNWoUKPftt9/qnnvuUb9+/dSnT59AmfHjx2vl\nypU6ffq0fv/736tz584xM8ePH9fw4cO1b98+/fjjj7rnnnt05ZVXBuo7ceKEunfvrnvuuUc9e/aM\n+efT09N1//33q3bt2pKkOnXq6E9/+lOgrlmzZunll19WQkKCBg0apA4dOkT982+//bZmzZp19uuM\njAytWrUqZs/Ro0c1bNgwHTp0SKdOndK9996rtm3bRs3k5ubq4Ycf1nfffafChQtr1KhRqlWrVtRM\n3p/vDz/8oAceeEA5OTmqWLGinnjiCRUpUiRqRpJef/11jRs3Tp9//rlKlCgRuOvBBx/U6dOnlZCQ\noCeeeEIVK1aMmlm1apXGjx+vhIQEFSlSRE888YTKlSsXs+uMJUuWqH///tqwYUPMzPDhw/XVV1+p\nTJkykqQ777zT9+edN3fq1CkNHz5c27ZtU4kSJTRhwgSVLl06Zm7QoEE6cOCAJOngwYNq0qSJxowZ\nEzXzxRdf6Omnn1ZCQoISExM1fvz4QF2bNm3SQw89pLi4OKWlpWnUqFFKSPj5oSnv87hhw4YxZ8Mv\n17lz55jz4dcVazb8chUrVgw0H5GOUdHmI29m4cKFgeYjb+7KK6+MOR95M7Nnz445G365smXLBpqP\nvLlatWpFnQ+/43W9evVizkek43y0+YjUFWs+/HJlypSJOh/RXoeizYZfbv78+VHnwy/Tpk2bmLPh\nl3vnnXdizodfLikpKep8+GWqVasW89hxxrmvza1atQp0/Mj7eh7ktcWvK8jx49xMjRo1Ah07/NYo\nRZ8Pv9znn38e6PhxbqZHjx6BXlvy5hYtWhTo+JE3l5qaGuj4cW6mcePGMefD7z1Y//79o85HpPdt\nseYjUles+fDLde/ePeqMRHtvGW0+/HJHjx6NOh9+meHDh8ecD7/cnj17Ys6HX65r165R58Mv07t3\n70DHj7zvtevWrRvo+PEz3gWWnp7u/c///I/neZ63ceNG7ze/+U2g3NGjR70+ffp4I0eO9KZMmRIo\ns3z5cq9///6e53ne/v37vfbt2wfKzZkzx3vxxRc9z/O8rKwsr3PnzoFynud5Tz/9tNezZ09vxowZ\ngf78Z5995g0cODDwf/+M/fv3e507d/aOHDni7dq1yxs5cqRTPj093Rs1alSgPztlyhTvySef9DzP\n83bu3Ol16dIlZubDDz/07r//fs/zPG/btm1nf+aR+P18hw8f7s2dO9fzPM976qmnvH/84x8xM++8\n84739NNPex06dPCys7MDdz3wwAPenDlzPM/zvKlTp3rjxo2LmRk4cKCXmZnpeZ7nTZw40Zs0aVKg\nLs/zvBMnTnh9+vTxrrjiikCZYcOGeQsXLoz2EPrmpk6d6o0ZM8bzPM+bNm2at2DBgsBrPGP48OHe\nmjVrYmZuuOEGb9OmTZ7ned6kSZO8yZMnB+oaMGCAt2jRIs/zPO+5557zZs2a9bOM3/M41mxEysWa\nD79MrNmIlAsyH5GOUdHmwy8TZD78crHmI9Yx1G82IuWCzIdfLtZ8+B2vg8yHXy7WfPhlgsyHXy7W\nfER6HYo2G5FysebDLxPk2BHrtTLSfPjlYs2HXybWbJzr3NfmIPORNxPktcUvF2Q+8maCHDv8cp4X\nez78ckGOH3kzQeYj0hrPiDQffrkgx4+8mSDz4fceLNZ8+GWCzIdfLsh8+OVizUik95ax5sMvF2s+\n/DJB5iPW+99I8+GXizUffpkg8+H3Xjvo8eNcF/wywuXLl+vqq6+WJNWqVUuHDh1SdnZ2zFyRIkX0\n0ksvKTk5OXBX8+bN9eyzz0qSSpUqpePHjysnJydmrlu3brrrrrskST/88IMqVaoUqG/Tpk3auHFj\nzDNM58Py5cvVqlUrJSUlKTk5OeK/DEXy/PPP65577gn0Z8uWLauDBw9Kkg4fPqyyZcvGzGzduvXs\nGctq1appx44dUR97v59venq6rrrqKknSlVdeqeXLl8fMXH311Ro8eLDi4uKcuh5++GF16dLl3/6+\n0TITJkxQamqqPM/Trl27VLly5UBdkvTCCy+od+/evv8aYpn1SLl//etfuu666yRJN99889nHM2jf\n5s2bdeTIkX87++yXOfdxO3TokO+c+OW2bdt29r/ftm1bLV269GcZv+dxrNmIlLvqqquizodfJtZs\nRMr95S9/iTkfkY5R0ebDelzzy8Waj2hdkWYjUq506dIx58Mvd+6xxG8+/I7XQebDLxfr+OGXCTIf\nfrlYx49Ir0PRZiNaLhq/TJBjR7SuaPPhl4t1/PDLxDp2nJH3tTnIfOTNBHlt8csFmY+8mSCvLX45\nKfZ8RMrFkjcTZD6idUWbD79ckNeXvJmg85FXkPnIK+h85BVkPvwEnZG8gszH+RB0PiKJNR95BZmP\nvILMh997bct8XPDN1t69e3/2IJQrV0579uyJmUtISFCxYsWcuuLj45WYmChJmj59utq1a6f4+PjA\n+VtuuUVDhgzRiBEjAv35cePGafjw4U5rlKSNGzdqwIAB+u1vfxv4YJCVlaUTJ05owIAB6t27d6Af\n9hlr165VlSpVfC9l8HPttddqx44d6tSpk/r06aNhw4bFzNSpU0effvqpcnJytHnzZm3fvv3sqWE/\nfj/f48ePnz0glC9f/t/mxC+TlJQUc21+ucTERMXHxysnJ0dvvPGGevToETMjSZ988om6du2qvXv3\nnj2wxMpt2bJF69ev1zXXXBN4fZI0depU3XbbbRo8eLD2798fKPf999/rk08+Ud++fTV48GDfA3m0\n59brr7/ue8muX2bEiBG699571aVLF61cuVI33HBDoFydOnW0ePFiST9d3rB3796f/f9+z+NYsxEp\nV7JkSd+/Z7RMrNmIlIuPj485H365zMzMqPMRqSvWfPjlYs1HtGNopNmIlBs5cmTM+fDL1a1bN+p8\nnHHu8TrIfPjlghw/8maCzIdfTop9/MibiXXsiNYVaz7yZoIcOyJ1SdHnwy8X5PiRNxPr2HFG3tfm\nIPORNxN0NvLmgsyH33uHILORNxd0Pvz6Ys1H3kzQ+Yj0vijWfOTNBZmPvJmg85H3PViQ+cibCTof\neXNBjx9+7xNjzUjeTND58OuKNR95M0HnI9L731jzkTcXZD7yZoLMh997bZfXlzP+4zfI8Dwv9I4F\nCxZo+vTpeuihh5xy06ZN06RJkzR06NCY63z33XfVpEkTpaamOnWkpaXpvvvu06RJkzRu3Dj98Y9/\n1MmTJwNlDx48qOeee05jx47Vgw8+GPixnD59esQXMj/vvfeeqlatqo8++kh///vfNXr06JiZ9u3b\nq2HDhrr11lv197//XTVr1vxFP+sLMSc5OTl64IEHdPnll6tVq1aBMu3atdMHH3ygmjVr6sUXXwyU\nefzxx/Xggw86re3666/XkCFD9Prrr6t+/fp67rnnAuU8z1ONGjU0ZcoU1a5dW5MnTw7cefLkSa1c\nuVKXX355oD8/ZswYPffcc5o/f76aNWumN954I1Bu2LBhmjdvnm677TZ5nhfxZx3peRxrNizP/7yZ\noLORNxd0Ps7NBZ2PczMu83FuLuh85P17BZ2Nc3Mu83FuLuh8RDpex5oPl+N8pEzQ+cibCzIf52Zc\njh3n5oLOx7mZ3NzcwMeOvH+voPNxbm706NGB5uPcTJDZiPXabMlEEikXbT4iZWLNhl8uyHz45WLN\nh18myLEj0t8t1nz45WIdP/wyQebD7z3YuVcKBM0Eed8WKRfr+BEpF21G/DKPPvpozPnwy1177bVR\n5yPS+mLNR6RcrPnwyz3yyCNR58MvM3To0ECvLdHeawd9zbjgN8hITk7+2e5x9+7dgc+wWCxZskQv\nvPCCXn755Zj/qn1GRkaGypcvrypVqqh+/frKycnR/v37Vb58+YiZRYsWafv27Vq0aJF27typIkWK\nqHLlymrdunXUrkqVKqlbt26SfrrcrkKFCtq1a1fMg3z58uV16aWXKiEhQdWqVVOJEiVirvGM9PR0\njRw5MuafO+PLL79UmzZtJEn16tXT7t27lZOTE/Ms4eDBg8/+76uvvjrQ2s6VmJioEydOqFixYtq1\na5fzZXWuHnzwQVWvXl333XdfoD//0UcfqVOnToqLi1OXLl00ceLEmJldu3Zp8+bNGjJkiKSf5r9P\nnz6aOnVq1Ny5B+COHTtq1KhRgdZYoUIFNW/eXJLUpk2bQGs844svvgh8Cl+SNmzYoGbNmkmSWrdu\nrffffz9QrkqVKmcPxEuWLNHu3bv/7c/kfR4HnQ3L898vE2Q28uaCzse5uWPHjgWaj7xdQecjby7I\nfPg9HkFmI28u6HzkzZUsWTLqfPgdr0uUKBFzPizH+UiZcePGRZ0Pv9y8efPUrVu3iPORN3P06FFt\n3Lgx5mz4ddWpU+fs38tvPvwyhQoVijkbkR6P9evXR50Pv1x6enrU+fDLFClSJOaxw++1Odbxw/p6\nHin37rvvRpwPv0zx4sV1zTXXRD125M0lJCSoUKFCMefDr2/06NGqX7++JP/58MsEOXZEejw8z4s6\nH365w4cPR52PSF2x5sPvPdi6deuizof1fVuk3MSJE6MeP/xyb7zxhvr16xdxRvJmihQpoqVLl8ac\nD7+utLS0s383v/nwy+zbty/mfER6PDIzM6POh1/u22+/jToffpnc3NyY8+H3Xjs+Pt75vekFP7N1\nxRVXaP78+ZKkr776SsnJyYFPv7o6cuSIxo8fr8mTJ5+9i0oQK1as0Kuvvirpp8sejx07FvP6z2ee\neUYzZszQW2+9pZtuukn33HNPzAOz9NNdTl555RVJ0p49e7Rv375A19i3adNGn332mXJzc3XgwIFA\na5R+erNfokQJp+t1q1evrjVr1kj66dKBM8MWzfr168/+C8onn3yiiy++WIUKuY1b69atz87Khx9+\nGPMOiL/ErFmzVLhwYQ0aNChwZuLEifrmm28kSWvWrFGNGjViZipVqqQFCxborbfe0ltvvaXk5OSY\nGy1JGjhwoLZv3y7pp83ymTvqxNKuXTstWbJE0k/PtyBrPGPdunWqV69e4D9foUIFbdy48Wy2evXq\ngXITJkzQokWLJEkzZ85Ux44df/b/+z2Pg8yG5fnvlwkyG365IPORNxdkPvy6gsyHXy7WfER6DGPN\nhl8uyHz45WLNh9/xOsh8WI7zfpmlS5fGnA+/3KRJk6LOR95Mbm5uoGOHX9dDDz0UdT78Mtdff33M\nY0ekxzDWfPjlateuHXU+/DJTpkyJOhtS5NfmaPNhfT33y+3duzfqfPhlYs2GX+6+++4LNB9+ff/8\n5z+jzodfpmPHjjHnI9LjGGs+/HKVKlWKOh9+mRUrVsScD7/3YD179ow6H9b3bX65FStWxDx++OXe\nfPPNqDOSN3Pq1CmtWbMm5nz4dY0dOzbqfPhlevXqFXM+Ij2OsebDL1euXLmo8+GXeeutt2LOh997\nbct70zjvQlyflceTTz6pFStWKC4uTg8//HCgN3QZGRkaN26cvv/+eyUkJKhSpUqaOHFi1DdRb775\npiZOnPizH/K4ceNUtWrVqF0nTpzQH//4R/3www86ceKE7rvvPt8fQiQTJ07URRddFOjW79nZ2Roy\nZIgOHz6sU6dO6b777lP79u0D9UybNk3Tp0+XJN19992BfgExIyNDzzzzjF5++eVAHdJPt34fMWKE\n9u3bp9OnT+v++++PeZldbm6uRowYoY0bN6po0aJ68sknVaVKlajryvvzffLJJzV8+HD9+OOPqlq1\nqh5//HEVLlw4aqZ169ZatmyZVq9erYYNG6pJkyZ64IEHYnbt27dPRYsWPbvxr1Wr1s/+9cYvM3To\nUD322GOKj49XsWLFNH78+H/7V/FYc9uxY0ctXLgwZqZPnz568cUXVbx4cSUmJurxxx8P1PXkk0/q\n0Ucf1Z49e5SYmKhx48apQoUKgdY4ceJENWvW7Oy/BsXKDB48WOPHj1fhwoVVunRpPfbYYypVqlTM\n3JAhQzRmzBh5nqfLLrvs3y5z8Hsejx07ViNHjow4G5FyLVu2VHp6esT58Mvs2LFDpUqVijgbkXKD\nBg3SU089FXU+Yh2j/ObDL9OzZ09NnTo16nxE6ho7dmzE+YiUeeWVVyLORqzHI9p8+OUGDhyo8ePH\nR5wPv+N1gwYNNGzYsKjz4ZfbsGFD1OOHX+bFF1/Ujz/+GHU+/HIVK1bUo48+GnE+Yr0O+c1GpFxi\nYqKeeOKJiPPhl2nVqpWGDRsW9dgRaY1jxoyJOh9+uTO3wo80H36ZtLQ0PfDAAxFnI68zr81t2rSJ\nOR95M7t27Yr52uKXe+utt2LOR95M7dq1o85GpNy57zkizYdfrmrVqlHnwy9zzTXXxJyPSGuMNR9+\nubS0tJivL3kzZ35G0ebD7z1Y/fr1o86HX+brr7+OOR9+uUmTJsWcD79cuXLlos5IrPeWkebDL1e0\naNGo8+GXadGiRcz5iLTGWPPhlytZsmTU+fDLpKamBjp+5H2v3bBhw8DHjzP+I5stAAAAACjo/uM3\nyAAAAACAgojNFgAAAACEgM0WAAAAAISAzRYAAAAAhIDNFgAAAACEgM0WAAAAAISAzRYAAAAAhOD/\nANuJb29oeIIHAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 1098x378 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "metadata": {
        "id": "OER2wALRvTkK",
        "colab_type": "code",
        "outputId": "d6af4fa4-c1b1-4619-af53-7f802ea5764c",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 679
        }
      },
      "cell_type": "code",
      "source": [
        "def plot_w(mrf,i,j):\n",
        "  \n",
        "  n = int(np.where((mrf[\"w_idx\"][:,0] == i)&(mrf[\"w_idx\"][:,1] == j))[0])\n",
        "  w = mrf[\"w\"][n]\n",
        "  \n",
        "  mx = np.max((w.max(),np.abs(w.min())))\n",
        "  plt.figure(figsize=(states/4,states/4))\n",
        "  plt.imshow(-w,cmap='bwr',vmin=-mx,vmax=mx)\n",
        "  plt.xticks(np.arange(0,states))\n",
        "  plt.yticks(np.arange(0,states))\n",
        "  plt.grid(False)\n",
        "  \n",
        "  ax = plt.gca()\n",
        "  al_a = list(alphabet)\n",
        "  ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x,y: al_a[x])) \n",
        "  ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x,y: al_a[x]))\n",
        "  plt.title(f\"W for positions {i} and {j}\")\n",
        "  plt.show()\n",
        "\n",
        "for n in range(2):\n",
        "  i = int(top.iloc[n][\"i\"])\n",
        "  j = int(top.iloc[n][\"j\"])\n",
        "  plot_w(mrf,i,j)"
      ],
      "execution_count": 14,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAAFLCAYAAABLIjB/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XtYVVX+P/D34SaSigmoOWKAo6jQ\nQVNEwNJM5aJIomYT8ijmlPcmFaFCM+xilpfRzJrHfFS8lpAaIJVaxuOFEQVvMzqJBxUTMhEVRa7r\n90ffzk8QOHtv98bLfr+exz/2cn/O+pzD4cPal7W2QQghQESkA1b3OwEiosbCgkdEusGCR0S6wYJH\nRLrBgkdEusGCR0S6wYL3gBg2bBhSUlLM2+Xl5fDx8UFqaqq5raysDE899RTy8vJqxN64cQPh4eEY\nPHgwrl692lgpN2j27NnYs2cPACAtLQ0lJSV3tWvpq6++wpAhQxAcHIxXXnkFBQUF5v87f/48hg8f\njnHjxmmeR7du3ZCfn39Xe25uLqKiohASEoKwsDB8//335v9LSkpCaGgoQkJCEB0dDZPJpHmeuiHo\ngbBgwQLx9ttvm7czMzOFn5+fmDNnjrlt//794rnnnrsr9tChQ+LZZ59tlDyVCAoKEpcuXWq0/o4e\nPSoCAwNFYWGhEOKPz3bGjBlCCCFyc3NFcHCwmDNnjhg7dqzmuXTt2lVcuHDhrvYhQ4aIrVu3CiGE\nOHXqlOjevbu4fv26OHPmjOjdu7coKCgQQgixceNG8dJLL2mep15whPeACAgIwIEDB8zbBw8exMiR\nI5GZmVmjLSAgoEbcr7/+ilmzZuHKlSsIDg5GUVERMjMzMXz4cAQHB2PUqFE4fvw4ACA5ORlTp07F\n2LFjsXDhwrty8PT0xLp16xAeHg5/f39s2rTJ/H/r1q1DaGgogoODMWnSJBQVFQEA/v3vf2P48OHm\nEcnOnTsBAFFRUdi+fTvefPNNmEwmREVFISsry9wOoME8p0+fjrfeegtBQUEIDQ3FL7/80mB/d2rV\nqhWWLFmC1q1bAwB69eqFM2fOAACaNGmCtWvXonv37g3+PKqrq/Huu+8iKCgIAwYMQExMDCoqKgAA\ncXFxWLZsGaKjo/Hcc88hOjoapaWlAIC9e/di0KBBCAkJwapVq+p87aqqKkyePBnh4eHmz93W1hb5\n+fnIzc2Fm5sb2rRpAwDo06eP+b2TCu53xaU/lJaWCm9vb/No4G9/+5s4evSoGDRokPmv/YsvvihS\nU1Pvij148KAYOHCgEEKIkpIS4efnJ7KysoQQQqSnp4vBgweLqqoqkZSUJLp37y5MJlOdOXTu3Fkk\nJCQIIf4YCXl7e4uioiKRnZ0tnn32WfH7778LIYRISEgQb731lhBCiIiICJGZmSmEEMJkMplHUmPG\njBHbtm0zv+6fI7w/2y3l6ePjI44fPy6EEGLevHnm0W99/TVk7ty5NUbKQgiRlJTU4AgvPT1dDB06\nVJSXl4vbt2+LkJAQ8/uJjY0VISEh4urVq6KiokIMGzZMbN++XVRWVorAwECRkZEhhBDiyy+/FJ07\nd65zhHennJwc4evrK27duiUKCgpEnz59xOnTp0V1dbX4/PPPxbRp0yy+R5KGI7wHhL29PXr27IkD\nBw6gtLQUubm58PLygq+vLw4ePIiSkhKcPHkS/v7+Db7OsWPH0LZtW/Ts2RMAEBQUhKtXr+LixYsA\nADc3N7i5udUbP2LECACAh4cH3N3dcezYMfz0008ICgqCk5MTAGDUqFHYt28fAMDJyQnbtm0zj0wW\nLVok6f1ayrNjx47w9vYG8Md5sEuXLinqb9u2bcjIyMC0adMk5fWnoKAgJCUlwdbWFk2aNMFTTz2F\nCxcumP+/X79+aNmyJWxsbNC5c2dcunQJeXl5KC8vR9++fQEAw4cPt9jPpUuXMHPmTMTHx6Np06Zo\n06YNZsyYgRdeeAG9e/fGhg0bMGvWLFm5U/1Y8B4gAQEBOHjwII4cOQIfHx9YW1ujd+/eyMzMxOHD\nh9G5c2c8/vjjDb5GUVERWrRoUaOtefPmuHLlCgDA0dGxwfg7/9/R0RHXr1+/6zVbtGhhfr0PPvgA\nTZs2RXR0NAYPHoz09HRJ79VSns2bNze3W1tbo6qqSnZ/GzZswIoVK7B27Vq4uLhIyuvO/GJjYxEU\nFITg4GDs3r0b4o5p53Xld+3aNTRr1szcbumzPnv2LKKiovDaa69h2LBhAID//Oc/WLlyJXbt2oVD\nhw5h5syZmDRpUo2+STkWvAdI3759cfjwYWRlZaF3794AAD8/P2RlZSErKwuBgYEWX8PJyQnFxcXm\nbSEErl27Zh6dWXLnVd7i4mI4OjrC2dm5xmsWFxfD2dkZAODs7Iw5c+bg559/xty5c/Hmm2/i5s2b\nmuUptb/k5GRs2LAB69evh6urq8V8aluyZAlsbGzw7bffIj09Hf369bMY4+joaL4aDcB8nrMuhYWF\nmDBhAmbOnIlRo0aZ2w8cOIAePXqgXbt2AIDQ0FCcOXPmgbn6/rBjwXuAdO3aFWVlZdi1axf8/PwA\nAG3btgXwx8nw2hcs6mI0GvH7778jOzsbAJCamoq2bduiffv2knL48zaY3NxcnDt3Dj4+Pujfvz9+\n+OEH8y/d5s2b0a9fP1RUVCAqKgq//fYbAMDLyws2Njawsqr5tbKxscH169fvOU+p/RUWFmLx4sVY\ntWqV+eS/XFeuXEHnzp1hZ2eHU6dOITs7G7du3WowpkOHDrC2tjZfaEpOTobBYKhz33feeQdjx45F\nSEhIjXZ3d3dkZ2ebP+u9e/fCxcXF4siepLG53wnQ/2cwGODv74+ff/4Z3bp1M7f7+voiJSXFfL6r\nIQ4ODli6dCnmz5+PW7duoVWrVli8eHG9v3i1tWrVCuHh4SgsLER8fDwcHR1hNBrx6quvIjIyEtXV\n1ejatSvmzZsHW1tbjBw50nw/m5WVlflc1J2Cg4Px0ksv4b333runPKX2t23bNty8eRPjx483t9nY\n2CAlJQWbNm3C2rVrUVJSgpKSEgQHB8NoNN511Xr8+PGIjY1FcnIyevXqhdjYWLz99tswGo0N5jd/\n/ny89dZbsLOzQ0REBBwcHO7ar7CwED/++CNMJlONK+GzZ8/GgAEDcPLkSbz00ksAgGbNmmHp0qWS\nf37UMIPgyQH6P56enti7d695VEn0qOEhLRHpBgseEekGD2mJSDc4wiMi3dD8Ku0dN6fL4no0xfJO\ntT31lKK+yto+KTumiZ2ygXFFpbKrbXfcsiZZq1aKuoJ1yTVlgXfcdCtVWaW1oq6aFJyTH/R/97bJ\nZqPw16SyUnbI4WO2yvpSqOfTjXiAV14uP8bOTllf9VzV5giPiHSDBY+IdIMFj4h0Q3bBS0lJgZeX\nV4PzBImIHkSKCp6rqyu+++47LfIhItKMrIJXXFyMY8eOIS4ursazFoiIHgayCl56ejr69++PZ555\nBnl5eSgsLNQqLyIi1ckqeCkpKRg6dCisra0RHByMtLQ0rfIiIlKd5DsqCwoKcPToUSxYsAAGgwG3\nb99G8+bNER0drWV+RESqkVzwUlJSEBkZibi4OAB/rFA7ePBgnD9/Hh06dNAsQSIitUg+pE1NTUVE\nRIR522Aw4IUXXuDFCyJ6aEge4X3zzTd3tU2ZMkXVZIiItMSZFkSkG9qvh/d/j92TTcFKE7C3V9SV\nqajhx+nVxd2tcVdL2b9ffky/gApFfcFK4d9BpXFKHDkiP+bpp9XPQ213PPVMjquVzS3vVIfHWzbe\nainXrsv/7ju2UJgfV0shIr1jwSMi3WDBIyLdkHyVNj8/H2FhYfD29gYAlJeXo3Pnzpg3bx6srZWt\nWktE1JhkjfDc3d2RmJiIxMREbNmyBRUVFfj222+1yo2ISFX3dEhrNBpx7pyCZwsQEd0HigteRUUF\ndu/eDS8vLzXzISLSjKzHMZlMJkRFRQEATp8+jQkTJmDgwIGaJEZEpDZZBe/Pc3gAMH36dLi7u2uS\nFBGRFhQf0sbExOCTTz5BaWmpmvkQEWlGccFzdXVFUFAQVq5cqWY+RESakXxI2759eyQnJ9domzFj\nhuoJERFpRfvFA27cUBan4HaXCk9vRV0pWaegqb2yj+2/p5QtHqBkXYRWrRR1BUcHhYsO2Mg6JXxv\nqqvlxzTm4gZKKVw8QNHnAQAtWiiLU6K8XH6MnZ2yvrh4ABHpHQseEekGCx4R6Ybkgpefn4+uXbvi\n1KlT5rbk5OS7LmQQET2oZI3w/vrXv2LRokVa5UJEpClZBc/LywsODg44cOCAVvkQEWlG9jm8N954\nA0uXLoXWd7MQEalNdsFzc3NDt27dkJaWpkU+RESaUXSVdsqUKfjXv/6FSiV37BIR3SeKCp6zszMG\nDhyIzZs3q50PEZFmFN+HN378eBQUFKiZCxGRpjiXFpxLWxvn0t5HnEtbE+fSEhEpo/2fZAcHZXHN\nm8sOOX9eWVcdy/8rP6hLF0V9uboqCkOzfd/JD3r2WWWdNcLXwkzJX30Ab8Q1kR2zZPFDcCvVrVvK\n4po1UzcPDVwpkf8zc2ql7s+MIzwi0g0WPCLSDRY8ItINySdr8vPzERYWBm9vbwghYG1tjYkTJ8Lf\n31/L/IiIVKP4MY3nz5/HxIkTsXjxYnRReAKfiKgxKT6k7dChAyZOnIiNGzeqmQ8RkWbu6Ryet7c3\nzpw5o1YuRESauqeCd/PmTVhbW6uVCxGRpu6p4J04cQJdu3ZVKxciIk0pLnjnz5/HmjVrMG7cOBXT\nISLSjqyrtCaTCVFRUSgvL0dVVRXmzp2Ldu3aaZUbEZGqJBe89u3bIzs7W8tciIg0pf3yUHz2BRE1\nNi4PRUR6x4JHRLrBgkdEuiHrKm1eXh4++OADFBUVobq6Gj169EBsbCzslC7DTETUiCSP8KqqqjBt\n2jRMmDABW7duRVJSEgBgxYoVmiVHRKQmySO8ffv2wcPDA7179wYAGAwGxMTEwOpheDAKERFkFLyz\nZ8/eNY3MXsmjtIiI7hPJwzODwYCqqiotcyEi0pTkgufh4YHjx4/XaCsvL8f//vc/1ZMiItKC5IIX\nGBiIixcvYs+ePQCA6upqfPzxx0hLS9MsOSIiNcmaWvbbb79h7ty5+O2332BnZ4eAgABMnTq14QsX\nnFpGRI2tnqllnEtLRI8ezqUlIr2TNdNCkfJyZXFFRbJDrtg9oagrJSk+0VbhyFXp56EkzsFBUVeX\ni5Qt2+/i3Iij+evX5ce0aKF+Hiq7Wlz3yMSSxx3KlHXYmLOkqqvlx6h8ny9HeESkGyx4RKQbLHhE\npBuSz+FduHAB77//Pi5fvgwhBHr16oWZM2eiSZMmWuZHRKQaSSO86upqTJs2DWPHjkVSUhKSk5PR\ntm1bzJ07V+v8iIhUI6ng7du3D08++ST8/f3NbdHR0cjJyUGRgqupRET3g6SCd/bsWXTr1q1Gm8Fg\nQKdOnWAymTRJjIhIbZIKnhCizpVShBDQeqIGEZFaJBU8d3d3nDhxokabEAJnzpyBh4eHJokREalN\nUsHr27cvcnNzsXfvXnPbmjVr0KNHD7Rq1Uqz5IiI1CR58YALFy4gNjYWJSUlEEKgR48eePvtty3f\nllKmcMoLp5bdexynltXEqWV3e1Snlt3r4gGurq7YuHEj5s2bh6ZNm2LevHm8B4+IHiqyy+fTTz8N\no9GIiIgI7Ny5U4uciIg0of16eLWWhZfKYHxKdoy4rXBY/5//yI/p3l1ZXzk5yuIWLpQfs2aNoq6u\n3lI2cre1VRSmyE8/yY/p319ZX80eU/YrcqNE/uFpc9xQ1Bfy85XFdemiLE6J4mL5MS1bKuuL6+ER\nkd6x4BGRbrDgEZFuyFrxOD8/H2FhYfD29q7Rvnz5crRUeqxNRNRIZC/x7u7ujsTERC1yISLSFA9p\niUg3WPCISDdkH9KaTCZERUWZt93d3ZGQkKBqUkREWuA5PCLSDR7SEpFu3PMhLQDExMTAaDSqlhQR\nkRZkFbz27dsjOztbq1yIiDSl+eIBly8ri3O5fUF2jGjvqqivs2flx3T0UPaxXbuubL0zx7yj8oNq\nPYdEMqVrkCmJU7JGGoCySvlr9jWxewgeR6DwoVi5xU6K4pR+jxVRsqaj0vX6uHgAEekdCx4R6QYL\nHhHphuSLFufOncOHH36IK1euAADatWuHd955hw/xIaKHhqQRXlVVFaZNm4YJEybg66+/xtdffw0v\nLy+8//77WudHRKQaSSO8ffv2oVOnTujVq5e5bcKECXwINxE9VCQVvLNnz8LT07NGm5XSWxeIiO4T\nSQXPysoKlZWV5u1JkyahpKQEBQUF2LFjB5o2bapZgkREapE0TOvUqROO3/H0sZUrVyIxMRFVVVWo\nVnjjKBFRY5NU8Pr06YOCggLs2bPH3Hby5EncvHkT1tbKnlJPRNTYJB3SGgwGrFq1CgkJCVixYgVs\nbW3h4OCAlStXwt7eXusciYhUIfk+PCcnJ/zzn//UMhciIk3xUisR6Ybmq6WA9+oRUWPjailEpHcs\neESkGyx4RKQbkgtefn4+IiIiarQtX74c69evVz0pIiItcIRHRLrBgkdEuiHrqWW1H9F48eJFjB8/\nXvWkiIi0IKvgubu7IzEx0by9fPly1RMiItIKD2mJSDdY8IhIN1jwiEg3OJeWiB49nEtLRHon6ypt\nYyq5WXeFbkizxxSOJpUsU9/YDzG6fl1+TIsW6ufxoLjjGSuS2TywX3ezy7/L/94DgIszj6Sk4AiP\niHSDBY+IdOOeCl5dCwoQET2oOMIjIt1gwSMi3WDBIyLdYMEjIt1gwSMi3WDBIyLdYMEjIt14YBcP\n4NSyWji1rCZOLauBU8tq4eIBRKR3D+wIr1H/gj8Eo4XC3+T/5W/T+hH+q/8Q/MzoPuIIj4j0jgWP\niHSDBY+IdEN2weMKKUT0sOIIj4h0gwWPiHSDBY+IdIMFj4h0gwWPiHSDBY+IdEPRXBuTyYSoqCjz\ndkxMDIxGo2pJERFpgXNpG7svhTiXtpaH4GdG9xHn0hKR3mk+wrt8WVmcy89J8oP69FHUV5nzX2TH\nNLFT9rFVVStb76x7d/kxx7PKFPWleK0/JSMoJWsRAsBXX8mPeeklZX01JiXrHgLKf2bNmimLU0BA\n/nffAIXliSM8ItI7Fjwi0g1ZBa/2wgG7du1CZGQkysvLVU+MiEhtii9bnT59GsuWLcOaNWtgZ2en\nZk5ERJpQdEhbVFSE2NhYLFmyBK1atVI7JyIiTcgueJWVlZg+fTpCQkLQsWNHLXIiItKE7IJnMpkQ\nEhKCpKQkFBQUaJETEZEmZBe8Tp06ITIyEjNmzMCsWbNQVVWlRV5ERKpTfFtKcHAwXF1dsWLFCjXz\nISLSzD3dhxcfH4+0tDRkZmaqlQ8RkWZk3ZbSvn17JCcnm7cfe+wxpKenq54UEZEWONOCiHTjwV0e\niohIKS4eQER6x4JHRLrBgkdEuiGp4OXn58PT0xM5OTk12keMGIG4uDhNEiMiUpvkEZ6rqytSUlLM\n2+fOncN1pauzEhHdB5ILno+PD/bv32+eSpaamorAwEDNEiMiUpvkgmdrawsfHx/zrIrdu3ejX79+\nmiVGRKQ2WTMtgoODkZKSAmdnZ7Rp0wYODg5a5UVEpDpZV2n9/f2RmZmJ1NRUBAUFaZUTEZEmZBU8\nOzs7+Pr6IikpCQMGDNAqJyIiTch+pkVwcDCKiorQvHlzLfIhItIM59IS0aOHc2mJSO9Y8IhIN1jw\niEg3WPCISDdY8IhIN2TdlpKfn4+wsDB4e3ub27p06YK3335b9cSIiNQm+z48d3d3JCYmapELEZGm\neEhLRLrBgkdEuiH7kNZkMiEqKsq8HRAQgEmTJqmaFBGRFngOj4h0g4e0RKQbLHhEpBtcLYWIHj1c\nLYWI9I4Fj4h0gwWPiHSDBY+IdEPx4gFCCJSXl+Pvf/87Bg0apFV+RESquacbj4uLizF8+HA888wz\nsLe3Vz05IiI13dMhbcuWLeHi4oLLly+rlQ8RkWbuqeDl5+ejuLgYTzzxhFr5EBFpRvHiAUIINGnS\nBB999BFsbGS/DBFRo+PiAUSkG7wthYh0gwWPiHSDiwcQ0aOHiwcQkd5pfnm18Le6K60lbeyuyo4p\nc3hcUV9NrCrkBym9Mn3rlrK46mrZIVfKmyvqyqlllaI4WDXi38/KSvkxD8HdBNeuK/t9cbQvU9ah\nnZ2yOCUUfIfV/k5xhEdEusGCR0S6wYJHRLoh66TGhg0bsH37dtjZ2eH27duYMWMGAgICtMqNiEhV\nkgtefn4+vvrqK2zduhW2trbIy8tDfHw8Cx4RPTQkH9KWlJSgrKwMFRV/XNF0c3PD+vXrNUuMiEht\nkgtely5dYDQa8fzzzyMuLg5paWmoVHJrABHRfSJ7pkVubi4yMjKwY8cOPPbYY1i3bh0M9dzVDACF\nhcoS4314tfA+vJp4H14NvA+vlnudaSGEQFlZGTp27Ihx48bh66+/RmFhIX799VdlCRERNTLJBW/r\n1q2YM2cO/hwQ3rhxA9XV1XByctIsOSIiNUke40dERODs2bMYNWoUHBwcUFlZifj4eD7LgogeGpqv\nlsJzeLXwHJ46eA6vBp7Dq4WrpRCR3j2w6+FdLZb/l+7xlgrfygOwigPJ9IiO8EglHOERkd6x4BGR\nbsga4y9YsAAnT57E5cuXUVpaig4dOsDR0RGffvqpVvkREalG0Tm85ORk/PLLL4iNjbW8M8/hkRZ4\nDo8awnN4RKR3LHhEpBsseESkGyx4RKQbLHhEpBsseESkG5xaBvC2lIcRb0uhhvC2FCLSuwd2hEdE\npBhHeESkdyx4RKQbLHhEpBsNFrywsDCcP3/evB0aGoq9e/eat6dMmYKMjAztsiMiUlGDBc/Pzw+H\nDh0CABQVFaG0tNS8DQBHjx5Fz549tc2QiEglFgteVlYWAODIkSMYNmwYcnJyAPzxQO727dvDwcFB\n+yyJiFTQYMHz9fXF4cOHAQBZWVkICAhAVVUVbt++jUOHDsHPz69RkiQiUkODBa9ly5ZwcHBAYWEh\njh49Ch8fHxiNRuTk5CArKwt9+vRprDyJiO6Zxau0fn5+yMjIgMFggL29PXr27Ins7GwcP34cPXr0\naIwciYhUIangbdmyBd27dwcA9OzZEz/99BNcXFxgb2+veYJERGqxWPB8fX1x8uRJ89VYJycnFBcX\n8/wdET10OJeWiB49nEtLRHrHgkdEusGCR0S6wYJHRLrBgkdEuiG54H322WdYvHixebu6uhrh4eE4\ndeqUJokREalNcsEbP348vvvuOxQWFgIAkpKS4OPjgy5dumiWHBGRmiQXPHt7e0yePBlLly5FaWkp\nVq9ejddff13L3IiIVCXrxmMhBEaPHg1XV1d4enri1VdflRJ0L/kREclXz43HsmdaHDhwADExMdiz\nZw/s7OwsB7DgEVFjU2umhaurK1q3bi2t2BERPUB4WwoR6QYLHhHpBldLIaJHD1dLISK9s9G8h+pq\nRWHZx6xlx3h4KOoKjnal8oOUrvZcWaks7vff5ce0bausr4fB9evyY1q0UD8PtSn9flgpHLsojVOg\nrLzuUVdDmtipe4TIER4R6QYLHhHphuSCN3r0aJw4caJG26JFi7B69WrVkyIi0oLkgjd06FDs3Lmz\nRtv333+PIUOGqJ4UEZEWJBe80NBQ/PDDD+btEydOoHXr1mjTpo0miRERqU1ywXNycoKrqyuOHTsG\nANi5cyfCwsI0S4yISG2yLloMHToUaWlpAIA9e/YgKChIk6SIiLQgq+ANGjQIP/74I44fPw43Nzc4\nOjpqlRcRkepkFbxmzZrB09MTX3zxBQ9nieihI/s+vLCwMOzbtw8DBgzQIh8iIs3Inlo2aNAgZGdn\na5ELEZGmONOCiHSDy0MR0aOHy0MRkd6x4BGRbrDgEZFusOARkW6w4BGRbrDgEZFusOARkW6w4BGR\nbih+atnGjRuxc+dOPP7441i2bJmaORERaYIzLYjo0cOZFkSkdyx4RKQbLHhEpBuKL1pIVs+xNBFR\nY+MIj4h0gwWPiHSDBY+IdIMFj4h0gwWPiHSDBY+IdIMFj4h0474UvJSUFHh5eaGoqEjS/vn5+ejR\noweioqIQFRWF0aNHY86cOaiqqrIY17VrV5w6dcrclpycjOTkZMn9jRkzBmPHjsWBAwcs5pmXl4dX\nX30VI0eOREREBObPn4/y8vIGYy5cuICJEydixIgRiIiIwAcffICysjKLfdX+TP78V1xcXG/MuXPn\nMHHiRIwaNQqjRo3C66+/bvFnkJ+fj4iIiBpty5cvx/r16y3mWF+8mvvXFbNr1y5ERkY2+Nnn5+fD\n09MTOTk5NdpHjBiBuLg4i/3V/uzff/99yTFjxozBiy++iB9++MHie9uwYQNefPFFjBkzBiNHjsT+\n/fstxixYsABRUVEIDg5Gv379EBUVhalTp9a7f1hYGM6fP2/eDg0Nxd69e83bU6ZMQUZGRp2xn332\nGRYvXmzerq6uRnh4eI3fudpGjx6NEydO1GhbtGgRVq9ebfG93TNxH7z22msiKChIbNy4UdL+Fy5c\nEMOHD6/RFhsbK7755huLcUOHDhUTJkwwtyUlJYmkpCRZ/Z07d06EhISI//73v/XGVFZWiqFDh4rM\nzEwhhBDV1dUiISFBLF68uN6YqqoqER4eLvbv329u+/LLL8Xs2bMbzK+uHC2prKwUYWFh4tChQ+a2\nL774QsyYMUN2P8uWLROJiYmS+pWbp9z9a8ecOnVKhIWFiStXrliMef7558X8+fPNbXl5eWLgwIEi\nNjZW0xyFEOLq1auif//+orS0tMGYYcOGifLyciGEECaTSURGRkruMykpSSxYsMDifvPnzxdbt24V\nQghx5coV0b9/f/Hxxx+b/z8wMFDcvHmzztjS0lIxePBgUVBQIIQQ4quvvhJz5sxpsL9169aJhQsX\n1mi78zW01OgjvOLiYhw7dgxxcXFITU1V/DpGoxHnzp2zuJ+XlxccHBwkjdDq06FDB0ycOBEbN26s\nd599+/bBw8MDvXv3BgAYDAbExMRgypQpDcY8+eST8Pf3N7dFR0cjJydH8uhXqn379qFTp07o1auX\nuW3ChAlYuHChqv3cT0VFRYjO6gbmAAAGGElEQVSNjcWSJUvQqlUri/v7+Phg//795iOF1NRUBAYG\nap0mAKBly5ZwcXHB5cuX692npKQEZWVlqKioAAC4ublJHlnL4efnh6ysLADAkSNHMGzYMPPINzc3\nF+3bt4eDg0Odsfb29pg8eTKWLl2K0tJSrF69Gq+//nqD/YWGhtYY3Z44cQKtW7dGmzZtVHpH9Wv0\ngpeeno7+/fvjmWeeQV5eHgoLC2W/RkVFBXbv3g0vLy9J+7/xxhtYunQpxD0sVeXt7Y0zZ87U+/9n\nz55F165da7TZ29vDzs6uwZhu3brVaDMYDOjUqRNMJpPiXOvry9PTs0ablZUVrK2tLcaaTKYah2/f\nfPONqrmpobKyEtOnT0dISAg6duwoKcbW1hY+Pj7IzMwEAOzevRv9+vXTMk2z/Px8FBcX44knnqh3\nny5dusBoNOL5559HXFwc0tLSUFlZqXouvr6+OHz4MAAgKysLAQEBqKqqwu3bt3Ho0CH4+fk1GD9s\n2DDk5uYiPj4ew4cPh5OTU4P7Ozk5wdXVFceOHQMA7Ny5E2FhYeq8GQu0n0tbS0pKCiZPngxra2sE\nBwcjLS0N0dHRFuP+/KUDgNOnT2PChAkYOHCgpD7d3NzQrVs3pKWlKc775s2bDRYHg8Fg8ZxibUKI\nOmOEEJKK852fCQC4u7sjISGhzn2trKxq/LJMmjQJJSUlKCgowI4dO9C0adN6+3F3d0diYqJ5e/ny\n5RZza2wmkwlxcXFYu3YtwsPD0bZtW0lxwcHBSElJgbOzM9q0aVPvSKau/u787AMCAjBp0iRJMUII\nNGnSBB999BFsbBr+FVy4cCFyc3ORkZGBVatWYdOmTVi3bh0MKs5Rb9myJRwcHFBYWIijR4/iH//4\nB4xGI3JycpCVlYURI0Y0GG8wGPDGG28gJiYGH374oaQ+hw4dirS0NBiNRuzZswebN2+WnO+9LD7c\nqAWvoKAAR48exYIFC2AwGHD79m00b95cUsG785du+vTpcHd3l9X3lClT8MorryAyMtLil6wuJ06c\nuGsEdycPDw9s2LChRlt5eTny8vLQuXPnOmPc3d2xZcuWGm1CCJw5cwYeHh4Wc6pdiBrSqVMnrFu3\nzry9cuVKAMCAAQNQXV0t6TUeZJ06dUJkZCScnJwwa9YsrF27VtLo1d/fHwkJCXBxcUFQUJDk/uR8\n9kpjhBAoLy9Hx44d0bFjR0RFRSEkJAS//vor/vKXv8jq2xI/Pz9kZGTAYDDA3t4ePXv2RHZ2No4f\nP4733nvPYryrqytat27d4BHNnQYNGoTPP/8cQ4YMgZubGxwdHSXn+vLLL+Pll1+WvP+dGvWQNiUl\nBZGRkdixYwe2b9+O9PR0XLt2rcYVIiliYmLwySefoLS0VHKMs7MzBg4cKOsvyZ/Onz+PNWvWYNy4\ncfXuExgYiIsXL2LPnj0A/rha9fHHHzc4quzbty9yc3NrXBFbs2YNevToIekclBx9+vRBQUGBOT8A\nOHnypMWR68MmODgYrq6uWLFihaT97ezs4Ovri6SkJAwYMEDj7OTZunUr5syZYx7t37hxA9XV1RYP\nGZXw8/PDli1b0L17dwBAz5498dNPP8HFxQX29vaq99esWTN4enriiy++aLTDWaCRR3ipqan46KOP\nzNsGgwEvvPACUlNTLR4O3MnV1RVBQUFYuXIlZsyYITlu/Pjx2LRpk6R9/zz8KC8vR1VVFebOnYt2\n7drVu7+VlRW+/PJLzJ07F59++ins7OwQEBDQ4O0A1tbWWLVqFWJjY7Fo0SIIIdCjRw+8++67snK8\nU0xMDIxG4137GgwGrFq1CgkJCVixYgVsbW3h4OCAlStXavKFvhe131d976k+8fHxGDFiBPz8/Cye\nfwL+KJJFRUVo3ry5ony1EhERgbNnz2LUqFFwcHBAZWUl4uPjNfl5+fr6YurUqZg4cSKAP86zFRcX\nY+jQoar39aewsDDMnj0bn3zyiWZ91Kb9My1IkiNHjmDBggXYvHkzrKx4PziRFvib9YB4+umnYTQa\nERERgZ07d97vdIgeSRzhEZFucIRHRLrBgkdEusGCR0S6wYJHRLrBgkdEuvH/AEtfP08X8PyLAAAA\nAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 378x378 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAAFLCAYAAABLIjB/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XtcVHX+P/DXcBMIwwQv+RUDXUSF\nQENFwNK8cAs08VYhj9R1Xc1Lm0pQqWtYZt7yoZm1D/NhmqYlmAZI62Uz1wsrihfc1U0EdEzURFSU\n68zn90fr/By5nXOcMwrn9Xw8/GM+nvd83jPMvOdzLp/P0QkhBIiINMDmUSdARGQtLHhEpBkseESk\nGSx4RKQZLHhEpBkseESkGSx4j4khQ4YgLS3N9LiyshIBAQFIT083tVVUVODZZ59FQUGBWezt27cx\ndOhQhIWF4caNG9ZKuV5vv/029u7dCwDIyMhAaWlpjXY13blzBzNnzkS3bt1q/N+uXbsQFhaGgQMH\nYtq0aabc1NCtWzfo9foa7QMGDMCgQYMQERFh+ndPSkoKoqKiEBkZiXHjxiE/P1+1/DRH0GNh4cKF\n4r333jM9zsrKEkFBQWLOnDmmtoMHD4oXX3yxRuyRI0fECy+8YJU8lQgPDxeXL1+2ap/R0dFi6dKl\nomvXrmbtFy5cEKGhoaKgoEAYjUbxwQcfiB07dqiWR9euXcXFixdrtL/44ou1tp87d0707t1bFBUV\nCSGE2LRpk3jllVdUy09rOMJ7TISEhODQoUOmx4cPH8aIESOQlZVl1hYSEmIW9+uvv2LWrFm4fv06\nIiIiUFxcjKysLAwbNgwREREYOXIkTp06BQBITU3F1KlT8frrr2PRokU1cvDx8cH69esxdOhQBAcH\n45tvvjH93/r16xEVFYWIiAhMnjwZxcXFAIB//etfGDZsmGlEsnPnTgBAfHw8tm/fjnfeeQf5+fmI\nj49Hdna2qR1AvXlOnz4d7777LsLDwxEVFYVffvml3v4elJycjFGjRtVo37FjB8LCwvDMM89Ap9Ph\nvffeQ0xMTI3tjEYj3n//fYSHh2PAgAFISEhAVVUVACApKQkrVqzAuHHj8OKLL2LcuHEoKysDAOzb\ntw+DBw9GZGQk1qxZU2tu9cnLy4OnpyfatGkDAOjTp4/ptZMFPOqKS78rKysTfn5+pl/9V199VZw4\ncUIMHjzY9Gs/atQokZ6eXiP28OHDYtCgQUIIIUpLS0VQUJDIzs4WQgiRmZkpwsLChMFgECkpKaJ7\n9+4iPz+/1hw6d+4skpOThRBC5OXlCT8/P1FcXCxycnLECy+8IH777TchhBDJycni3XffFUIIERsb\nK7KysoQQQuTn54sZM2YIIYQYM2aM+P77703Pe2+Ed6+9oTwDAgLEqVOnhBBCzJs3zzT6rau/2ly8\neLHGCG/atGliwYIFYuzYsSIsLEzMmTNH3L17t0ZsZmamiI6OFpWVlaK8vFxERkaaXk9iYqKIjIwU\nN27cEFVVVWLIkCFi+/btorq6WoSGhor9+/cLIYT48ssvRefOnesc4U2bNk1ER0eL2NhYsXv3biGE\nEEVFRaJPnz7i7Nmzwmg0is8//1xMmzatztdI8nCE95hwdHREYGAgDh06hLKyMuTl5cHX1xe9evXC\n4cOHUVpaitOnTyM4OLje5zl58iTatm2LwMBAAEB4eDhu3LiBS5cuAQA8PT3h6elZZ/zw4cMBAB07\ndoSXlxdOnjyJn376CeHh4XBzcwMAjBw5EgcOHAAAuLm54fvvvzeNTJYuXSrp9TaUZ6dOneDn5wfg\n9+Ngly9ffqj+7rl16xYOHjyIJUuWYNu2bbh48SI+//zzGtuFh4cjJSUF9vb2aNasGZ599llcvHjR\n9P/9+vVDixYtYGdnh86dO+Py5csoKChAZWUl+vbtCwAYNmxYnXlERUUhLi4OP/zwA9555x0kJCSg\nsLAQbdq0wYwZM/Dyyy+jd+/e2LhxI2bNmiXrNVLdWPAeIyEhITh8+DCOHTuGgIAA2Nraonfv3sjK\nysLRo0fRuXNnPPXUU/U+R3FxMZ588kmztubNm+P69esAAFdX13rj7/9/V1dX3Lp1q8ZzPvnkk6bn\nW7BgAZycnDBu3DiEhYUhMzNT0mttKM/mzZub2m1tbWEwGB6qv/v7GDhwINzc3ODs7IxXX33VVLwf\nzC8xMRHh4eGIiIjAnj17IO6bdl5bfjdv3oSLi4upvb73etasWQgKCgIA9OzZE71798Y///lP/Pvf\n/8bq1auxe/duHDlyBDNnzsTkyZPN+iblWPAeI3379sXRo0eRnZ2N3r17AwCCgoKQnZ2N7OxshIaG\nNvgcbm5uKCkpMT0WQuDmzZum0VlD7j/LW1JSAldXV7i7u5s9Z0lJCdzd3QEA7u7umDNnDn7++WfM\nnTsX77zzDu7cuaNankr7u6ddu3ZmZ2VtbGxga2tbY7tPPvkEdnZ2+OGHH5CZmYl+/fo1+Nyurq5m\nz33vOOeDKisraxyXMxgMsLe3x6FDh9CjRw+0a9cOwO8jwXPnzj02Z98bOxa8x0jXrl1RUVGB3bt3\nm37927ZtC+D3g+EPnrCojb+/P3777Tfk5OQAANLT09G2bVu0b99eUg73LoPJy8tDYWEhAgIC0L9/\nf+zatcv0pdu8eTP69euHqqoqxMfH4+rVqwAAX19f2NnZwcbG/GNlZ2eHW7duPXSeUvurT2RkJDIy\nMlBUVASDwYCtW7fWepjg+vXr6Ny5MxwcHHDmzBnk5OTg7t279T53hw4dYGtrazrRlJqaCp1OV2O7\nsrIyjB492vTaz549i2PHjiE4OBheXl7Iyckxvdf79u1Dq1atGhzZkzR2jzoB+v90Oh2Cg4Px888/\nm10/1qtXL6SlpZmOd9XH2dkZy5cvx/z583H37l20bNkSy5Ytq/WLV5uWLVti6NChuHLlCmbPng1X\nV1f4+/tj4sSJiIuLg9FoRNeuXTFv3jzY29tjxIgRGDt2LIDfR0uzZ8+Gk5OT2XNGRETglVdewQcf\nfPBQeUrt7/Tp05g5cyaqq6thMBhM17hlZmaie/fumDp1Kl577TXY2dkhMDAQEydOrNHX+PHjkZiY\niNTUVPTs2ROJiYl477334O/vX29+8+fPx7vvvgsHBwfExsbC2dm5xnaurq5Yvnw5/vrXv6KiogJO\nTk5YvHgxPDw84OHhgdOnT+OVV14BALi4uGD58uWS/35UP53gwQH6Hx8fH+zbt880qiRqarhLS0Sa\nwYJHRJrBXVoi0gyO8IhIM1Q/S6t0GuAD16RK4uiorC/XY/+QH9S/v7LOfvtNWVwDl0TU6n/Xyskm\n4zKPh45zcFDU1e1S+Wctm9uVKepL8QdLgdzTys7Gensr66+ZgxV38MrLZYforzs1vFEt6rq6iSM8\nItIMFjwi0gwWPCLSDNkFLy0tDb6+vnXOEyQielwpKngeHh748ccf1ciHiEg1sgpeSUkJTp48iaSk\nJLN7LRARNQayCl5mZib69++P559/HgUFBbhy5YpaeRERWZysgpeWlobo6GjY2toiIiICGRkZauVF\nRGRxki88LioqwokTJ7Bw4ULodDqUl5ejefPmGDdunJr5ERFZjOSCl5aWhri4OCQlJQH4fYXasLAw\nXLhwAR06dFAtQSIiS5G8S5ueno7Y2FjTY51Oh5dffpknL4io0ZA8wtu2bVuNtilTplg0GSIiNXGm\nBRFphvrr4f3vLvOyKViNRDxZ/y0I66JkkQ9hVPi2PXAzG6lyzst/bRLv21NDK8fbygIVrIahaFkc\nAPjuO9khVaPHKOrK3k7h37qoSH7Mfbd5lMPg3LzhjWpha2PF1VKUrBR0313gZPHyqrWZIzwi0gwW\nPCLSDBY8ItIMyWdp9Xo9YmJi4OfnB+D3u6d37twZ8+bNq/XO7UREjxtZIzwvLy9s2LABGzZswJYt\nW1BVVYUffvhBrdyIiCzqoXZp/f39UVhYaKlciIhUpbjgVVVVYc+ePfD19bVkPkREqpF117L8/HzE\nx8cDAM6ePYsJEyZg0KBBqiRGRGRpsgrevWN4ADB9+nR41XFxHxHR40jxLm1CQgKWLFmCsjKF9/ok\nIrIyxQXPw8MD4eHhWL16tSXzISJSjeRd2vbt2yM1NdWsbcaMGRZPiIhILeovHmAwKAoTNvIvZtZB\n4UuprJQf4+CgrC+FSu/oZMe4PGHFieFWVlYu//1wcmy670eTVV2tLM7evtZmTi0jIs1gwSMizWDB\nIyLNkFzw9Ho9unbtijNnzpjaUlNTa5zIICJ6XMka4f3hD3/A0qVL1cqFiEhVsgqer68vnJ2dcejQ\nIbXyISJSjexjeG+99RaWL18Ota9mISKyNNkFz9PTE926dUNGRoYa+RARqUbRWdopU6bgb3/7G6qV\nXhRIRPQIKCp47u7uGDRoEDZv3mzpfIiIVKP4Orzx48ejSMl9N4mIHhHOpQU4l7YR4lxajeBcWiIi\nZWSteKyIktETAJ2jo+yY26Xyf/UBoLmNgl8RhSO8m7eU5Wg0Kgprsqw8wFZGyR9N6R/aTv2v8iPB\nER4RkTIseESkGSx4RKQZknf89Xo9YmJi4OfnByEEbG1tMWnSJAQHB6uZHxGRxSi+TeOFCxcwadIk\nLFu2DF26dFElOSIiS1K8S9uhQwdMmjQJmzZtsmQ+RESqeahjeH5+fjh37pylciEiUtVDFbw7d+7A\n1lb+jAgiokfhoQpebm4uunbtaqlciIhUpbjgXbhwAevWrcPYsWMtmA4RkXpknaXNz89HfHw8Kisr\nYTAYMHfuXLRr106t3IiILEpywWvfvj1ycnLUzIWISFXqLw/Fe1+YUbLMEwC4OClYZsuGE2lIo3S1\nf8/4jSAizWDBIyLNYMEjIs2QdZa2oKAACxYsQHFxMYxGI3r06IHExEQ4NIrVGIlI6ySP8AwGA6ZN\nm4YJEyZg69atSElJAQCsWrVKteSIiCxJ8gjvwIED6NixI3r37g0A0Ol0SEhIgA3PBBJRIyG54J0/\nf77GNDJHBfedICJ6VCQPz3Q6HQwKb7lIRPQ4kFzwOnbsiFOnTpm1VVZW4r///a/FkyIiUoPkghca\nGopLly5h7969AACj0YjFixcjIyNDteSIiCxJ1tSyq1evYu7cubh69SocHBwQEhKCqVOn1n/iglPL\nzHBqGZEV1DG1jHNprYwFj8gKOJeWiLRO1kwLJfLOKxvRdOpQZeFM6pGVJT8mNFRRVyUlisLgcl0v\nP6htW0V9lRmbKYq7cEF+jMIU4eIiP8bWqOwzVVZtryjOCWWyY6rsnBT1ZW+j8AoKa+4FVFdbJwYA\nnGp/HznCIyLNYMEjIs1gwSMizZB8DO/ixYv48MMPce3aNQgh0LNnT8ycORPNmik73kNEZG2SRnhG\noxHTpk3D66+/jpSUFKSmpqJt27aYO3eu2vkREVmMpIJ34MABPPPMMwgODja1jRs3DsePH0dxcbFq\nyRERWZKkgnf+/Hl069bNrE2n08Hb2xv5+fmqJEZEZGmSCp4QotaVUoQQUHuiBhGRpUgqeF5eXsjN\nzTVrE0Lg3Llz6NixoyqJERFZmqSC17dvX+Tl5WHfvn2mtnXr1qFHjx5o2bKlaskREVmSpMtSbG1t\nsWbNGiQmJmLp0qUQQqBHjx54//331c6PiMhiJF947OHhgU2bNmHevHlwcnLCvHnzeA0eETUqsmda\nPPfcc/D390dsbCx27typRk5ERKpQfT28vDxlcS1ayI9xs7uprLPUVPkxY8cq6+vcOWVxa9bIDqn6\n4GNFXSleoMKmQn6QnbIFe6JibGXHZKQr/KgbjYrChI38HHW3lH2GS21dFcW5PGG9qywKL8hfOUnp\nfcLatKm9nXNpiUgzWPCISDNY8IhIM2QdQNHr9YiJiYGfn59Z+8qVK9FCyUE3IiIrkn3E2MvLCxs2\nbFAjFyIiVXGXlog0gwWPiDRD9i5tfn4+4uPjTY+9vLyQnJxs0aSIiNTAY3hEpBncpSUizXjoXVoA\nSEhIgL+/v8WSIiJSg6yC1759e+Tk5KiVCxGRqlRfPODaNWVxrdytuHR8ZaX8GAcHZX0pnIiO8nL5\nMc7OyvoiekxUVctfcAAA7O1rb+cxPCLSDBY8ItIMFjwi0gzJJy0KCwvx0Ucf4fr16wCAdu3a4a9/\n/Stv4kNEjYakEZ7BYMC0adMwYcIEfPfdd/juu+/g6+uLDz/8UO38iIgsRtII78CBA/D29kbPnj1N\nbRMmTOBNuImoUZFU8M6fPw8fHx+zNhsbHv4josZFUsGzsbFB9X13dpk8eTJKS0tRVFSEHTt2wMnJ\nSbUEiYgsRdIwzdvbG6dOnTI9Xr16NTZs2ACDwQCj0gtpiYisTFLB69OnD4qKirB3715T2+nTp3Hn\nzh3Y2sq/FR0R0aMgeWrZ9evXkZycDL1eD3t7ezg7O2Pq1Kl47rnn6o3j1LIHcGoZkWSWnlrGubQA\nCx7RY4pzaYmIFJK9Hp5cjeLqlbt35ccoHeEpVGH3hOyYZmjC10lac1RuTUr3ABrFF02+igplcRzh\nEZHmseARkWaw4BGRZkgueHq9HrGxsWZtK1euxNdff23xpIiI1MARHhFpBgseEWmGrMtSHrxF46VL\nlzB+/HiLJ0VEpAZZBc/LywsbNmwwPV65cqXFEyIiUgt3aYlIM1jwiEgzWPCISDNUXy3lfzc5k82t\npRXngZaUyI9p0UJZXwrnSlZUy193sJkD59Ka4VzaRqf0jrLVUlxcam9vmu8SEVEtVB/hlZUpi3Ny\nbKKjk0awHp7SNciUpNjcRdnfWUmO9nZN9DPVlCn5UAFAHffZ4QiPiDSDBY+INOOhCl5tCwoQET2u\nOMIjIs1gwSMizWDBIyLNYMEjIs1gwSMizWDBIyLNYMEjIs3g1DJr49QyM5xaRvXi1DIiImVkLfGu\nhJODQWGk9WqxwSh/tGBr03RHC0pHQvZ1LMmjBnsbJZ8r6/6+Xy6S/7l6um3T/VwpYuElvTjCIyLN\nYMEjIs1gwSMizZBd8LhCChE1VhzhEZFmsOARkWaw4BGRZrDgEZFmsOARkWaw4BGRZiiaWpafn4/4\n+HjT44SEBPj7+1ssKSIiNai+WgoMCufS2jTRubSNYLWURkHJ+2jFzxTAubQWofT7YmtbazN3aYlI\nM9Qf4eXmKovz8pIfY6ds8ZdtGc1kxwx7WeHbdvWqorDBY9rIjsnIUNQV7MtvKwtUMoJydFTW188/\ny48JCVHWl9IVO27dkh9TXKysr7ZtlcUpff+VuHtXfozS9fDc3Gpt5giPiDSDBY+INENWwXtw4YDd\nu3cjLi4OlZWVFk+MiMjSFK94fPbsWaxYsQLr1q2Dg4VXJSUiUoOiXdri4mIkJibik08+QcuWLS2d\nExGRKmQXvOrqakyfPh2RkZHo1KmTGjkREalCdsHLz89HZGQkUlJSUFRUpEZORESqkF3wvL29ERcX\nhxkzZmDWrFkwKJ1JQURkZYovS4mIiICHhwdWrVplyXyIiFTzUNfhzZ49GxkZGcjKyrJUPkREqpF1\nWUr79u2RmppqevzEE08gMzPT4kkREamBMy2ISDPUXzxA5acnIqpBV/vSXBzhEZFmsOARkWaw4BGR\nZkgqeHq9Hj4+Pjh+/LhZ+/Dhw5GUlKRKYkREliZ5hOfh4YG0tDTT48LCQtxSsqIrEdEjIrngBQQE\n4ODBg6apZOnp6QgNDVUtMSIiS5Nc8Ozt7REQEGCaVbFnzx7069dPtcSIiCxN1kyLiIgIpKWlwd3d\nHW3atIFzU74NIBE1ObLO0gYHByMrKwvp6ekIDw9XKyciIlXIKngODg7o1asXUlJSMGDAALVyIiJS\nhex7WkRERKC4uBjNmzdXIx8iItVwLi0RNT2cS0tEWseCR0SawYJHRJrBgkdEmsGCR0SaIeuyFL1e\nj5iYGPj5+ZnaunTpgvfee8/iiRERWZrs6/C8vLywYcMGNXIhIlIVd2mJSDNY8IhIM2Tv0ubn5yM+\nPt70OCQkBJMnT7ZoUkREauAxPCLSDO7SEpFmsOARkWZwtRQianq4WgoRaZ3skxayGY3K4mysWIvL\ny+XHODpaPo963LxV+y9WfVyfbMKjayWfK2t+phS6XSr/7wwAzV2a6N+6slJZXLNmtTY//p8AIiIL\nYcEjIs1QvHiAEAKVlZX405/+hMGDB6uVHxGRxTzUhcclJSUYNmwYnn/+eTha+ZgWEZFcD7VL26JF\nC7Rq1QrXrl2zVD5ERKp5qIKn1+tRUlKCp59+2lL5EBGpRvHiAUIINGvWDB9//DHs7NS/uoWI6GFx\n8QAi0gxelkJEmsGCR0Saof7iAQaDsjhOLTPDqWUP4NQyM5xa9gBOLSMirVP99GpFta2iuGYO8n+x\nSu8o+3V0sbde3b9cpCzHp5+8Iz/IaOWLwa05grp7V36Mi4vl86iPglFoc0eFi21YYR2Qh3XlqvzP\nfhtHBXtfAEd4REQseESkGSx4RKQZsnb8N27ciO3bt8PBwQHl5eWYMWMGQkJC1MqNiMiiJBc8vV6P\nb7/9Flu3boW9vT0KCgowe/ZsFjwiajQk79KWlpaioqICVVVVAABPT098/fXXqiVGRGRpkgtely5d\n4O/vj4EDByIpKQkZGRmorq5WMzciIouSPdMiLy8P+/fvx44dO/DEE09g/fr10NVxSzQAqKhQlph1\nr8NTkKSDg6K+rHodnrUXZbXmdXilpfJjGsF1eIpvetUIVixSdh3eTWWdubrW2iz5EyqEQEVFBTp1\n6oSxY8fiu+++w5UrV/Drr78qS4iIyMokF7ytW7dizpw5uDcgvH37NoxGI9zc3FRLjojIkiSPg2Nj\nY3H+/HmMHDkSzs7OqK6uxuzZs3kvCyJqNCQXPFtbWyQmJqqZCxGRqjjTgog0Q/VTO81KrigLdHeX\nHeLipKwrnPy3/Jju3RV1pWSRDwBA5W/yY9q1U9TVtRJ7RXGtWipY+1Dpmd1bt+THODsr60thjr+c\nl79SkIuLstWFnm79+K872cZdfo55BbWfbW1IpzrCOMIjIs1gwSMizZC1S7tw4UKcPn0a165dQ1lZ\nGTp06ABXV1d8+umnauVHRGQxsgpeUlISACA1NRW//PILz9oSUaPCXVoi0gwWPCLSDBY8ItIMFjwi\n0gwWPCLSDBY8ItIMRVPLYmNjLZ0HEZHqOMIjIs2QvcS7bCo/PRFRDXXcdoIjPCLSDBY8ItIMFjwi\n0ox6C15MTAwuXLhgehwVFYV9+/aZHk+ZMgX79+9XLzsiIguqt+AFBQXhyJEjAIDi4mKUlZWZHgPA\niRMnEBgYqG6GREQW0mDBy87OBgAcO3YMQ4YMwfHjxwH8fkPu9u3bw1npstlERFZWb8Hr1asXjh49\nCgDIzs5GSEgIDAYDysvLceTIEQQFBVklSSIiS6i34LVo0QLOzs64cuUKTpw4gYCAAPj7++P48ePI\nzs5Gnz59rJUnEdFDa/AsbVBQEPbv3w+dTgdHR0cEBgYiJycHp06dQo8ePayRIxGRRUgqeFu2bEH3\n/92WMDAwED/99BNatWoFR0dH1RMkIrKUBgter169cPr0adPZWDc3N5SUlPD4HRE1OpxLS0RND+fS\nEpHWseARkWaw4BGRZrDgEZFmsOARkWZILnifffYZli1bZnpsNBoxdOhQnDlzRpXEiIgsTXLBGz9+\nPH788UdcuXIFAJCSkoKAgAB06dJFteSIiCxJcsFzdHTEG2+8geXLl6OsrAxr167Fm2++qWZuREQW\nJevCYyEERo8eDQ8PD/j4+GDixIlSgh4mPyIi+eq48Fj2TItDhw4hISEBe/fuhYODQ8MBLHhEZG2W\nmmnh4eGB1q1bSyt2RESPEV6WQkSawYJHRJrB1VKIqOnhailEpHV2qvdgNCqLs7FiLS4vlx+jdLVn\npe+Hkhyb8h3lFLyPl6/aKurq6bZW3EtpDN8XKyq9U/tIrSEuLrW3N813iYioFix4RKQZkgve6NGj\nkZuba9a2dOlSrF271uJJERGpQXLBi46Oxs6dO83a/v73v+Oll16yeFJERGqQXPCioqKwa9cu0+Pc\n3Fy0bt0abdq0USUxIiJLk1zw3Nzc4OHhgZMnTwIAdu7ciZiYGNUSIyKyNFknLaKjo5GRkQEA2Lt3\nL8LDw1VJiohIDbIK3uDBg/GPf/wDp06dgqenJ1xdXdXKi4jI4mQVPBcXF/j4+OCLL77g7iwRNTqy\nr8OLiYnBgQMHMGDAADXyISJSjfqLBxgMyuI4tcwcp5aZ49Qyc5xaZoZTy4hI87g8FBE1PVweioi0\njgWPiDSDBY+INIMFj4g0gwWPiDSDBY+INIMFj4g0gwWPiDRD8V3LNm3ahJ07d+Kpp57CihUrLJkT\nEZEqONOCiJoezrQgIq1jwSMizWDBIyLNUHzSQrI69qWJiKyNIzwi0gwWPCLSDBY8ItIMFjwi0gwW\nPCLSDBY8ItIMFjwi0oxHUvDS0tLg6+uL4uJiSdvr9Xr06NED8fHxiI+Px+jRozFnzhwYGrjnrV6v\nR9euXXHmzBlTW2pqKlJTUyX3N2bMGLz++us4dOhQg3kWFBRg4sSJGDFiBGJjYzF//nxUVlbWG3Px\n4kVMmjQJw4cPR2xsLBYsWICKiooG+3rwPbn3r6SkpM6YwsJCTJo0CSNHjsTIkSPx5ptvNvg30Ov1\niI2NNWtbuXIlvv766wZzrCvektvXFrN7927ExcXV+97r9Xr4+Pjg+PHjZu3Dhw9HUlJSg/09+N5/\n+OGHkmPGjBmDUaNGYdeuXQ2+to0bN2LUqFEYM2YMRowYgYMHDzYYs3DhQsTHxyMiIgL9+vVDfHw8\npk6dWuf2MTExuHDhgulxVFQU9u3bZ3o8ZcoU7N+/v9bYzz77DMuWLTM9NhqNGDp0qNl37kGjR49G\nbm6uWdvSpUuxdu3aBl/bQxOPwJ///GcRHh4uNm3aJGn7ixcvimHDhpm1JSYmim3btjUYFx0dLSZM\nmGBqS0lJESkpKbL6KywsFJGRkeI///lPnTHV1dUiOjpaZGVlCSGEMBqNIjk5WSxbtqzOGIPBIIYO\nHSoOHjxoavvyyy/F22+/XW9+teXYkOrqahETEyOOHDliavviiy/EjBkzZPezYsUKsWHDBkn9ys1T\n7vYPxpw5c0bExMSI69evNxiAPICOAAAHJElEQVQzcOBAMX/+fFNbQUGBGDRokEhMTFQ1RyGEuHHj\nhujfv78oKyurN2bIkCGisrJSCCFEfn6+iIuLk9xnSkqKWLhwYYPbzZ8/X2zdulUIIcT169dF//79\nxeLFi03/HxoaKu7cuVNrbFlZmQgLCxNFRUVCCCG+/fZbMWfOnHr7W79+vVi0aJFZ2/3PoSarj/BK\nSkpw8uRJJCUlIT09XfHz+Pv7o7CwsMHtfH194ezsLGmEVpcOHTpg0qRJ2LRpU53bHDhwAB07dkTv\n3r0BADqdDgkJCZgyZUq9Mc888wyCg4NNbePGjcPx48clj36lOnDgALy9vdGzZ09T24QJE7Bo0SKL\n9vMoFRcXIzExEZ988glatmzZ4PYBAQE4ePCgaU8hPT0doaGhaqcJAGjRogVatWqFa9eu1blNaWkp\nKioqUFVVBQDw9PSUPLKWIygoCNnZ2QCAY8eOYciQIaaRb15eHtq3bw9nZ+daYx0dHfHGG29g+fLl\nKCsrw9q1a/Hmm2/W219UVJTZ6DY3NxetW7dGmzZtLPSK6mb1gpeZmYn+/fvj+eefR0FBAa5cuSL7\nOaqqqrBnzx74+vpK2v6tt97C8uXLIR5iqSo/Pz+cO3euzv8/f/48unbtatbm6OgIBweHemO6detm\n1qbT6eDt7Y38/HzFudbVl4+Pj1mbjY0NbG1tG4zNz883233btm2bRXOzhOrqakyfPh2RkZHo1KmT\npBh7e3sEBAQgKysLALBnzx7069dPzTRN9Ho9SkpK8PTTT9e5TZcuXeDv74+BAwciKSkJGRkZqK6u\ntnguvXr1wtGjRwEA2dnZCAkJgcFgQHl5OY4cOYKgoKB644cMGYK8vDzMnj0bw4YNg5ubW73bu7m5\nwcPDAydPngQA7Ny5EzExMZZ5MQ1Qfy7tA9LS0vDGG2/A1tYWERERyMjIwLhx4xqMu/elA4CzZ89i\nwoQJGDRokKQ+PT090a1bN2RkZCjO+86dO/UWB51O1+AxxQcJIWqNEUJIKs73vycA4OXlheTk5Fq3\ntbGxMfuyTJ48GaWlpSgqKsKOHTvg5ORUZz9eXl7YsGGD6fHKlSsbzM3a8vPzkZSUhK+++gpDhw5F\n27ZtJcVFREQgLS0N7u7uaNOmTZ0jmdr6u/+9DwkJweTJkyXFCCHQrFkzfPzxx7Czq/8ruGjRIuTl\n5WH//v1Ys2YNvvnmG6xfvx46C85Rb9GiBZydnXHlyhWcOHECf/nLX+Dv74/jx48jOzsbw4cPrzde\np9PhrbfeQkJCAj766CNJfUZHRyMjIwP+/v7Yu3cvNm/eLDnfh1l82KoFr6ioCCdOnMDChQuh0+lQ\nXl6O5s2bSyp493/ppk+fDi8vL1l9T5kyBX/84x8RFxfX4IesNrm5uTVGcPfr2LEjNm7caNZWWVmJ\ngoICdO7cudYYLy8vbNmyxaxNCIFz586hY8eODeb0YCGqj7e3N9avX296vHr1agDAgAEDYDQaJT3H\n48zb2xtxcXFwc3PDrFmz8NVXX0kavQYHByM5ORmtWrVCeHi45P7kvPdKY4QQqKysRKdOndCpUyfE\nx8cjMjISv/76K/7v//5PVt8NCQoKwv79+6HT6eDo6IjAwEDk5OTg1KlT+OCDDxqM9/DwQOvWrevd\no7nf4MGD8fnnn+Oll16Cp6cnXF1dJef62muv4bXXXpO8/f2sukublpaGuLg47NixA9u3b0dmZiZu\n3rxpdoZIioSEBCxZsgRlZWWSY9zd3TFo0CBZvyT3XLhwAevWrcPYsWPr3CY0NBSXLl3C3r17Afx+\ntmrx4sX1jir79u2LvLw8szNi69atQ48ePSQdg5KjT58+KCoqMuUHAKdPn25w5NrYREREwMPDA6tW\nrZK0vYODA3r16oWUlBQMGDBA5ezk2bp1K+bMmWMa7d++fRtGo7HBXUYlgoKCsGXLFnTv3h0AEBgY\niJ9++gmtWrWCo6OjxftzcXGBj48PvvjiC6vtzgJWHuGlp6fj448/Nj3W6XR4+eWXkZ6e3uDuwP08\nPDwQHh6O1atXY8aMGZLjxo8fj2+++UbStvd2PyorK2EwGDB37ly0a9euzu1tbGzw5ZdfYu7cufj0\n00/h4OCAkJCQei8HsLW1xZo1a5CYmIilS5dCCIEePXrg/fffl5Xj/RISEuDv719jW51OhzVr1iA5\nORmrVq2Cvb09nJ2dsXr1alU+0A/jwddV12uqy+zZszF8+HAEBQU1ePwJ+L1IFhcXo3nz5oryVUts\nbCzOnz+PkSNHwtnZGdXV1Zg9e7Yqf69evXph6tSpmDRpEoDfj7OVlJQgOjra4n3dExMTg7fffhtL\nlixRrY8HqX9PC5Lk2LFjWLhwITZv3gwbG14PTqQGfrMeE8899xz8/f0RGxuLnTt3Pup0iJokjvCI\nSDM4wiMizWDBIyLNYMEjIs1gwSMizWDBIyLN+H8Yx1PTyzZr3AAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 378x378 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "metadata": {
        "id": "xv2anXSUX9kr",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## validating output\n",
        "(comparing to known output from GREMLIN_cpp version)"
      ]
    },
    {
      "metadata": {
        "id": "RBIvsIQ9X74s",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "!wget -q -nc http://files.ipd.uw.edu/krypton/4FAZA.out"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "M7fK03GFX_3q",
        "colab_type": "code",
        "outputId": "250b8a81-836f-4691-aa60-a830dc600782",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 347
        }
      },
      "cell_type": "code",
      "source": [
        "pd_test = pd.read_table(\"4FAZA.out\",sep=\" \")\n",
        "plt.scatter(pd_test[\"apc\"],mtx[\"apc\"])\n",
        "plt.show()"
      ],
      "execution_count": 16,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAFKCAYAAAAnj5dkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VPW9//H3ZJIMCROywIysFi4Q\nkFAwCCoGQTFxa70PHz+rBFzAulx7uXXlZzG0BrUs1x9Yb+2mrdW6p/aX28ujPx+iVmxVgogCCl4I\ncAvKmhmykCHJZDIzvz8wkWiWmZOZzJmZ1/MfMznzPfn02zDvnHO+iyUYDAYFAABMIyXWBQAAgM4I\nZwAATIZwBgDAZAhnAABMhnAGAMBkCGcAAEwmNdYFtHO5GqN6/tzcTNXVNUX1ZyQK+io09FNo6KfQ\n0E+hSaR+cjiyuj2WNFfOqanWWJcQN+ir0NBPoaGfQkM/hSZZ+ilpwhkAgHhBOAMAYDKEMwAAJkM4\nAwBgMoQzAAAmQzgDAGAyhDMAACZDOAMAYDKEMwAAJkM4AwBgMoQzAABd8Pr8qqlrktfn7/efbZqN\nLwAAMAN/IKCKt/dqa7VLtSe8yhtkU2G+Q/PmjpM1pX+uaQlnAABOU/H2Xr215WDH6+MnvB2vFxTn\n90sN3NYGAOBLXp9fW6tdXR7bWu3ut1vchDMAAF9q8HhVe8Lb5bG6xhY1eLo+FmmEMwAAX8q225Q3\nyNblsdysAcq2d30s0ghnAAC+ZEuzqjDf0eWxwvwhsqVZ+6UOBoQBAHCaeXPHSTr1jLmusUW5WQNU\nmD+k4/v9gXAGACQVr8+vBo9X2XZbl1fC1pQULSjO1zVzxvb4vmginAEASSHc+cu2NKucuZkxqJRw\nBgAkCTPMXw4VA8IAAAnPLPOXQ0U4AwASnlnmL4eKcAYAJDyzzF8OFeEMAEh4Zpm/HCoGhAEA4lpv\nU6PamWH+cqgIZwBAXAp3apQZ5i+HinAGAMQlo1OjYjl/OVQ8cwYAxJ14mxoVrj5dOa9cuVLbt2+X\nxWJRWVmZpkyZ0nHsxRdf1Lp165SSkqLJkydr2bJlfS4WAAAptKlRZr867onhK+fNmzfrwIEDqqio\n0IoVK7RixYqOYx6PR08//bRefPFFvfzyy9q3b5+2bdsWkYIBAIi3qVHhMhzOVVVVKi4uliSNHTtW\nDQ0N8ng8kqS0tDSlpaWpqalJbW1tam5uVnZ2dmQqBgAkvXibGhUuw7e13W63CgoKOl7n5eXJ5XLJ\nbrfLZrNp8eLFKi4uls1m03e+8x2NGTMmIgUDACDF19SocEVstHYwGOz42uPx6Mknn9Trr78uu92u\nhQsXateuXZo4cWK37XNzM5WaGt2/dByOrKieP5HQV6Ghn0JDP4WGfvpKS2ub6k54lTvIpgHpnaPq\n9H66a/45Pb43Xhn+X+F0OuV2uzte19TUyOE4dYth3759GjVqlPLy8iRJ06dP144dO3oM57q6JqOl\nhMThyJLL1RjVn5Eo6KvQ0E+hoZ9CQz+d0tvc5e76KVVSY0Oz4qkHe/pjzPAz56KiIq1fv16StHPn\nTjmdTtntdknSiBEjtG/fPrW0tEiSduzYodGjRxv9UQCAJNE+d/n4Ca+C+mrucsXbe2NdWr8yfOU8\nbdo0FRQUqLS0VBaLReXl5aqsrFRWVpZKSkp0yy236KabbpLValVhYaGmT58eyboBAAmmt7nL18wZ\n288VxU6fbs4vWbKk0+vTb1uXlpaqtLS0L6cHACSRUOYuj+znmmKFFcIAAKaQ6HOXw0E4AwBMIdHn\nLocjMcacAwASQiLPXQ4H4QwAMI142tYxmghnAIDpxMO2jtHEM2cAAEyGcAYAwGQIZwAATIZwBgDA\nZAhnAABMhnAGAMBkCGcAAEyGcAYAwGQIZwAATIZwBgD0yOvzq6auSV6fP9alJA2W7wQAdMkfCOil\nN6u1dY9b9Z5WDR5kU2G+Q/PmjpM1hWu7aCKcAQDf4A8E9PCzW/RFjafje8dPePXWloOSpAXF+bEq\nLSnwpw8A4BteemtPp2A+3dZqN7e4o4xwBgB04vX5ta3a3e3x2hMtavB4+7Gi5EM4AwA6afB4Vd9D\n+Gbb05Vtt/VjRcmHcAYAdJJttylvUPfhWzh+iGxp1n6sKPkQzgCATmxpVhXmO7o8Nspp14ISBoNF\nG6O1AQDfMG/uOEmnBn/VNrYoZ6BNZ+cP0YLi8Uyj6geEMwAkCa/PrwaPV9l2W6+3pa0pKVpQnK9r\n5owNuQ0ih3AGgATXl8VEbGlWOXMz+6lStCOcASCB+QMBPfTMhzroOtnxPRYTMT8eHABAgmrytukn\nv9vcKZhPx2Ii5sWVMwAkGH8goIq39+rd7Yfk9QW7fV/7YiLctjYfwhkAEkzF23s7blv3hMVEzIvb\n2gCQQLw+v7ZWu0J6L4uJmBfhDABx7vT9lhs8XtWe6H3d65GOgSwmYmLc1gaAONX+bHlrtUu1J7zK\nG2TTlHFDlJuVrtrG1m7bjXQMVPnNM1hMxMQIZwCIU19/tnz8hFcbPj6kUU57l+GcnpaioslDtaAk\nn2A2OcIZAOJQT8+Wm1p8urhwuD7ZV6u6xhbl2G2a+K1cLSgZr0xbWj9XCiMIZwCIQz09W65r9Oqy\nc8/UdXPHs/RmnOK+BgDEoZ62dczNGtARyM7cTII5DhHOABCHetrWsTCfKVLxjtvaABCnTt/Wsa6x\nRblZA1SYP6Tj+4hfhDMAmEyoWzuyrWPiIpwBwCS6mrccytaObOuYeAhnADCJruYts7VjcjI8IGzl\nypWaN2+eSktL9cknn3Q6duTIEc2fP1/f+9739OCDD/a5SABIVO1LbzZ4vN3OW2Zrx+Rj6Mp58+bN\nOnDggCoqKrRv3z6VlZWpoqKi4/jq1av1/e9/XyUlJXrooYd0+PBhDR8+PGJFA0C8++Yt7AE63u28\nZbZ2TDaGrpyrqqpUXFwsSRo7dqwaGhrk8XgkSYFAQB999JHmzp0rSSovLyeYAeBr2m9hHz/hVVDS\n8RMt3b63fd4ykoehK2e3262CgoKO13l5eXK5XLLb7aqtrdXAgQO1atUq7dy5U9OnT9d9993X6zlz\nczOVmhrdUYYOR1ZUz59I6KvQ0E+hoZ9OaWltU90JrzIHpOqTfcdDblc0dbhGDs+JYmXxJRl+nyIy\nICwYDHb6+tixY7rppps0YsQI3X777XrnnXd00UUX9XiOurqmSJTSLYcjSy5XY1R/RqKgr0JDP4WG\nfvrmLexse7rqPd3vGpVjT9eJk60d85avmnlm0vdhu0T6ferpjwxD4ex0OuV2uzte19TUyOE4tVJN\nbm6uhg8frjPPPFOSNHPmTO3Zs6fXcAaARPX1Udg9BfPgQQP04KLpava2MW85iRl65lxUVKT169dL\nknbu3Cmn0ym73S5JSk1N1ahRo7R///6O42PGjIlMtQAQZ3raPaorhflDlJWZzprYSc7QlfO0adNU\nUFCg0tJSWSwWlZeXq7KyUllZWSopKVFZWZmWLl2qYDCo/Pz8jsFhAJBseto9SpJy7TY1nPRqSE6G\npowdzNKbkNSHZ85Llizp9HrixIkdX3/rW9/Syy+/bLwqAEgQ7btHdTVN6vRb2GNHD1ZjQ3MMKoQZ\nsSsVAERRb7tHtd/CHpDOgo34Cr8NABBl7B6FcBHOANAHoewgxe5RCBfhDAAGGNlBit2jECrCGQAM\nYAcpRBMDwgAgTD3NXWYHKUQC4QwAYepp7nL7DlJAXxDOABCm9rnLXWEHKUQC4QwAYept7jIjsdFX\nDAgDAAOYu4xoIpwBwADmLiOaCGcA6APmLiMaeOYMAIDJEM4AAJgM4QwAgMkQzgCShtfnV01dEyt4\nwfQYEAYg4Z2+ScXxE17l2NNVOH6IFpTkd7tJBRBL/FYCSHjtm1Qc/3LJzXpPqzZsPayHn90ifyAQ\n4+qAbyKcASS0njap+KLGo5fe2tPPFQG9I5wBJLQGj7fjirkr29hFCiZEOANIaNl2m3Ls6d0erz/p\nZRcpmA7hDCCh2dKsKhw/pNvjeewiBRMinAEkvAUl+RrltHd5jF2kYEaEM4CEZ01J0YOLpuviaSOU\na7fJYpEGDxqg4ukj2UUKpsQ8ZwBJwZqSohsvnaDrLh7HLlIwPa6cAcSdvqz01b6LFMEMM+PKGUDc\nOH2lr9oTXuUNsqkw36F5c8ex0hcSCuEMIG60r/TV7vgJb8frBcX5sSoLiDj+1AQQF3pa6WsrC4kg\nwRDOAOJCg8er2m5W+qprbGEhESQUwhmAKX190Fe23aa8QV0vFpLLQiJIMDxzBmAqTd42vbB+t3Z9\nXqd6T6sGnzboqzDf0emZczsWEkGiIZwBmII/ENDLf92jv209JP9puziePuirfcGQrdVu1TW2KDdr\ngArzh7CQCBIO4QzAFF56s1obth7u9vjWareumTNWC4rzdc2csSwkgoTGM2cAMeUPBPT8G7v1t23d\nB7Mk1Z74atAXC4kg0XHlDCBmvD6/Xli/W+/vONrre7Pt6Qz6QtIgnAH0u/aVvj7eXaPaxtaQ2hSO\nZ9AXkgfhDKDfeH1+ueqa9NqmA9r0WU3I7UY57VpQwgpgSB6EM4Coax+JvfHTI2ppDfTe4EvpqRZd\n8O1hur4kn7WzkVQIZwBR1dLapmdf2xXSc+XTnT/pDC28YiK3spGUDIfzypUrtX37dlksFpWVlWnK\nlCnfeM/atWu1bds2Pf/8830qEkD8aX+uvH2vW676lpDbDWanKcBYOG/evFkHDhxQRUWF9u3bp7Ky\nMlVUVHR6z969e/Xhhx8qLS0tIoUCiC9f30EqFBdMHqobL5vA1TKSnqE/S6uqqlRcXCxJGjt2rBoa\nGuTxeDq9Z/Xq1brnnnv6XiGAuOL1+XXQ5dHHu0Mf8JWXZVPx9JG6+UpuYwOSwStnt9utgoKCjtd5\neXlyuVyy2+2SpMrKSp177rkaMWJEZKoEYHr+QEAvvFGtrdUunWjyhdyuaPJQ3cDVMtBJRAaEBYPB\njq/r6+tVWVmpZ555RseOHQv5HLm5mUpNje4/TocjK6rnTyT0VWjop1MaPV7dvvoteZrbQm6TYUvV\nJTNG6dZ/niyrlWfLEr9PoUqGfjIUzk6nU263u+N1TU2NHA6HJGnTpk2qra3V9ddfr9bWVn3++eda\nuXKlysrKejxnXV2TkVJC5nBkyeVqjOrPSBT0VWjop68Gfb390UEFgr2/X5LOm3SGvjPzW3LkZMiW\nZlVt7cnoFhkn+H0KTSL1U09/ZBgK56KiIj3xxBMqLS3Vzp075XQ6O25pX3755br88sslSQcPHtQD\nDzzQazADiE+hDvqySMob9NUOUozCBnpmKJynTZumgoIClZaWymKxqLy8XJWVlcrKylJJSUmkawRg\nQl6fX1urXb2+L8du073XTZGDjSqAkBl+5rxkyZJOrydOnPiN94wcOZI5zkCCavB4VXvC2+v7zps8\nVCOdif+MEIgk7i0BMCTbblPeoJ53ibJnpOoH/+ubCxQB6BnhDKBHxxuatfHTIzre0Nzp+7Y0qwrz\nHd22G+EYqDWLL2AkNmAAa2sD6FJzq08/+nVVp+lR9oxU/fsPZioj/dTKf/PmjpMkba12q66xRTlZ\nNo0emqUbLp2gHPZeBgwjnAF8w/GGZv34dx/I6+u8g5SnuU0/+nWVfn7XbEmSNSVFC4rzdc2csWrw\neJVttzHoC4gAwhlAh66ulr/O09ym4w3NGpyd0fE9W5pVztzM/igRSAo8DAIgr8+vmrom3f+rnoO5\n3e7P6/uhKiB5ceUMJDF/IKCX3tqjbdVu1Xl6nxbVbsKZOVGsCgDhDCQpfyCgh575UAdd4S2fac9I\n7XRLG0DkEc5AEvIHAlr+zIc6ZCCY//0HM6NUFYB2hDOQhF56szqsYE5PS9GKW8/jihnoJ4QzkGS8\nPr+27nH3/sYvfX1uM4DoI5yBBOf1+TvNQW7weFXvae2xzcyCM1QwOk8TzszhahmIAcIZSFD+QEAv\nvVmtrXvcqve0avAgmwrzHbr6wn/S4EE2He9m04qRjoH6/nfOYltHIIYIZyABNXnb9PCzm1VT19Lx\nveMnvB17LxfmO7rch3mEY6DKb55BMAMxRjgDCcQfCOiFN3bpb9uOdvuerdVuPXTLjI6va0+0KNue\nrsLxQ7SgJJ9gBkyAcAYSRKjTo2pPtMjT5GNNbMDECGcgATR5ffrpH7boaG1zr+/Ntqcr+8sdo1gT\nGzAnwhmIY/Uer154Y7d2/M9xtbYFQ2pTOH4IV8mAyRHOQBxqbWvTiuc+1hc1nrDaDR+SqQUl+VGq\nCkCkEM5AnPH6/Fr++806dtpI7FAMtFn10PfPZcAXEAcIZyBO+AMBvfLXPXrvkyPy+gJhtR04IFVr\n/+0CghmIE4QzEAcam1r1zGu7tG1v6MtuSlKKRZo1ZZhuvGwCwQzEEcIZMDFPs1cP/2GL3PWh77Us\nSempFhXmO3XjZfnKtLEmNhBvCGfAhE4tJlKtv207HHbb8yc5dcNlE5Vp4583EK/41wuYjD8Q0MPP\nbgl7JLYkzTl7qBZePikKVQHoT4QzYCKNTa363WufhR3MtrQUXTh1uObNHRelygD0J8IZMIHWtjY9\n8uxHOuw+qdCWEjnFlpqiaRMcuv5Sni0DiYRwBmLM0+zVvU9sVFsgnFiWzp90hhZeMZHVvoAERDgD\nMeIPBFTx9l79dcvBsK6WJWnutOGaX8wOUkCiIpyBGPD6/Hp+/W5t3NH91o5dcWTb9JNF02XPsEWp\nMgBmQDgD/cgfCOilN6v1cbVLDSd9YbW9cOpQ3XwFI7GBZEA4A/2kyevTij98pCO1TWG1GzzIpsJ8\nByOxgSRCOANR1uRt0wvrd+njarda20JfEzvVatGDi2bIkZPBoC8gyRDOQJT4AwE9v3633vvkiMIc\niK2hgzO0/OYZSk/lnyiQjPiXD0RBa1ub7vvF+zrZ4g+77YVTztDNVxZEoSoA8YJwBiLMHwjo3l9s\nVFOYwZyXZdO0CTxbBkA4AxFV7/Fq5XNb1NTSFnKbkY6B+sHVk5U3aADPlgFIIpyBiGhsbtWq5z/W\n0TBGYqenpqhoyjAtKB7PYiIAOiGcgT7wBwJ68c3demfrkbDaDc3L1I8XTmdbRwBd4pMBMOhkc6se\n/N3msOctD3dk6qGbz+VqGUC3DIfzypUrtX37dlksFpWVlWnKlCkdxzZt2qTHHntMKSkpGjNmjFas\nWKEUPoiQINpX+dq446i8vtDnLUvS7LOH6sZLJxLMAHpkKJw3b96sAwcOqKKiQvv27VNZWZkqKio6\njj/44IN67rnnNHToUN1555169913NWfOnIgVDcSKPxDQ8mc+1CHXybDanZGboZ8smsFtbAAhMfRJ\nUVVVpeLiYknS2LFj1dDQII/HI7vdLkmqrKzs+DovL091dXURKheInYM1jVpTsV0nTraG3CbFIs05\ne7gWlLCDFIDQGQpnt9utgoKvFknIy8uTy+XqCOT2/9bU1Oj999/XXXfdFYFSgdiobWzWj35dJX94\nd7CVabNq9R3ns4MUgLBF5B5bMPjNtQmPHz+uO+64Q+Xl5crNze31HLm5mUpNje4cT4cjK6rnTyT0\nleT3B/TLP27Vm1sOhtXOYpGKZ4zS4u+dLauVq2WJ36dQ0U+hSYZ+MhTOTqdTbre743VNTY0cDkfH\na4/Ho9tuu0133323Zs2aFdI56+rCG/EaLocjSy5XY1R/RqKgr6TGZp8e+E2VmryhLyaSnpqiaflD\ndMNlE5VpS1VtbXjPpRMVv0+hoZ9Ck0j91NMfGYbCuaioSE888YRKS0u1c+dOOZ3OjlvZkrR69Wot\nXLhQs2fPNnJ6ICa8Pr+O1jbpLxv3a+tul8K5iz00N0M/XjRdmba0qNUHIHkYCudp06apoKBApaWl\nslgsKi8vV2VlpbKysjRr1iz9+c9/1oEDB/SnP/1JkvTd735X8+bNi2jhQKT4AwG98tc9eu+TI2FP\njZJOLb9ZfvMMBnwBiBjDz5yXLFnS6fXEiRM7vt6xY4fxioB+9tz6XXp3+1FDbS+eNoLlNwFEHJMu\nkbRa29r04O82q6a+xVD7Z5cVKxDuEG4ACAF/7iMpHahp1B1r/m4omDPSU/SbJbM1OG9gFCoDAK6c\nkWSaW31a8ov31dwa/hVvikUq+vZQ3XQ5y28CiC7CGUnD09KqOx9/L+x2tjSL7is9W6Ocg9hvGUC/\nIJyR8Lw+v/5xuEGPvrwt7LbDhmSofNEMpafyTwVA/+ETBwnLHwjo5beq9fbHhw21nzV1qL5/xaQI\nVwUAvSOckbB+81879NFud+9v/JoUSbOnDdf1xfmRLwoAQkA4I+H0ZdDXffOmatzIHJ4tA4gpwhkJ\n5dO9Lv3sT58aavvLey9URjrLbwKIPcIZCeFo3UmVPfmBobYD0ix6/K4LGfQFwDT4NEJca/K26ck/\nf6pP/1FnqH3Rt4dq0RXMWwZgLoQz4pI/ENCv//NTfbznuKH2Z48frFu/O4ldpACYEuGMuNPQ5NU9\nP3/fcPtHbjtXIwbbe38jAMQI4Yy44mpo0o9+vclQ28FZaVrxLzN5tgzA9PiUQlzwBwL61X9+qq0G\nb2OvvP1cDc3jahlAfCCcYXr/OHJCj/xhi+H2P7uzSNmZtghWBADRRTjDtBpOtuj+X22Uz2+s/ZDs\nNP30Nm5jA4g/fGrBdPyBgJ7403Z98j/GpkeNHWbX4mumKMc+IMKVAUD/IJxhKq4TzfrRr6oMtz9/\nklO3//PkCFYEAP2PcIZpbN3j0hP/19jSm5I0e+ow3XjZhAhWBACxQTgj5upPenXfE+8raLD9LVdM\n0PRJQ9msAkDCIJwRM/5AQP/+3BbtPeoxfA42qwCQiAhnxMR/7z+u//PKdsPtzz/LoVuuKmBNbAAJ\niXBGv2pta9NPfrdZrvoWQ+3TrNKj/3qBsgcyEhtA4iKcEVVen18NHq+y7Ta1+Nr6tCb2w7fM0EhH\nVgSrAwBzIpwRFf5AQBVv79XWapdqT3gND/ZqxypfAJIJ4YyoqHh7r97acrDP55k8epDuLZ0egYoA\nIH4Qzog4r8+vj3fX9Pk8Rd92aNEVBRGoCADiC+GMiPIHAvr1f25XbWOr4XOMG5qpu+efo0wbU6QA\nJCfCGRHT2Nyqu/7jvT6dY+W/nKehuQMjVBEAxCfCGX3mDwT05LpPtWWXsb2WJenWKyfqginDI1gV\nAMQvwhl90tzq0+LH3u3TOX6zZDbbOgLAafhEhGGfHajVmpe3GW6fa7dq1R1FBDMAfA2fighba1ub\n7nzs72oNGD/HAzcUavzI3MgVBQAJhHBGWD7aVaNf/nmH4fajnWlatqiINbEBoAeEM0LiaWnVnY/3\nbST2Azeeo/EjsiNUEQAkLsIZvXr348N65o1dhtsPzbHokdvncLUMACEinNGtSFwt/8ddFyorg8VE\nACAchDO6tPMfx7W2wvh+y+eelafbrprC1TIAGEA4oxN/IKB5S/5LTX3YRmr1Dy6QM5v9lgHAKMIZ\nHd7/5Iiefu2/DbfPybRo9b9eyLxlAOgjw5+iK1eu1Pbt22WxWFRWVqYpU6Z0HNu4caMee+wxWa1W\nzZ49W4sXL45IsYiO1rY23f3zv6vF+F4Vuv/6aZo4KidyRQFAEjMUzps3b9aBAwdUUVGhffv2qays\nTBUVFR3Hf/rTn+rpp5/WGWecoRtuuEGXXXaZxo0bF7GiETmbdhzVU3/5rE/nYPlNAIgsQ5+oVVVV\nKi4uliSNHTtWDQ0N8ng8stvt+uKLL5Sdna1hw4ZJkubMmaOqqirC2WRqGpq19NdVfTrHkEFWrbrj\nQgZ9AUCEGQpnt9utgoKCjtd5eXlyuVyy2+1yuVzKy8vrdOyLL77o9Zy5uZlKTbUaKSdkDkdWVM8f\nD/z+gB783d/1SXVDn85TfN4o/ds1Z8tqTe5g5ncqNPRTaOin0CRDP0XkXmQw2IehvV+qq2uKQCXd\ncziy5HI1RvVnmF0kdpC69YqzdM4kp2xpVtXWnoxQZfGJ36nQ0E+hoZ9Ck0j91NMfGYbC2el0yu12\nd7yuqamRw+Ho8tixY8fkdDqN/BhEUNWnh/Xb/2d8lS9J+u39F3ELGwD6gaFP2qKiIq1fv16StHPn\nTjmdTtntdknSyJEj5fF4dPDgQbW1tWnDhg0qKiqKXMUIS/3JVn1/9dt9CuYZE7IJZgDoR4aunKdN\nm6aCggKVlpbKYrGovLxclZWVysrKUklJiZYvX6777rtPknTllVdqzJgxES0avfMHAlr7YpV2HfIa\nPscgq7TyztnKtDESGwD6kyUYiQfGERDtZwiJ9JyiN7v2H9ejrxhfelOS/vf8s3XWt/J6f2MSS6bf\nqb6gn0JDP4Umkfop4s+cYU4NTV7d8/P3+3SOGeNz9OAdcxLmlx8A4hHhnCCaW319DuYHbjpH44ez\n3zIAxBrhnADe+uCAXtqwr0/n+OW9Fyojna0dAcAMCOc4dqS2Scue2tSncywpPVuTRvNsGQDMhHCO\nQ61tbbr/l3/XiWbj57j2wm/piqKxkSsKABAxhHOcOVp3UmVPftCncyy/9VydOcQeoYoAAJFGOMcJ\nfyCgn72yTZ99Xm/4HIMzpNU/ZDERADA7wjkORGKK1D3XTtW3xw6OUEUAgGginE2u6pND+u1ruw23\nnzY+W/92zTkRrAgAEG2Es0lF4mr5vtKzVcBIbACIO4SzCf3j6Ak98uwWw+2HDErVqjtm8WwZAOIU\n4WwirW1tuvvxv6ulzfg5fnLzdI05Y1DkigIA9DvC2STe2nRAL71jfJWvASnS4/fOVnoq/5cCQLzj\nkzzGDh336Ce/3dynczx867kaybxlAEgYhHOM+AMBPbluh7bschs+x8KSCZpzzogIVgUAMAPCOQb8\ngYDu/Nk7avYZP8fP7ixSdqYtckUBAEyDcO5nH+w4oif/8t+G248fatP9N81kJDYAJDDCuZ94Wlp1\n1+PvKdiHczASGwCSA+HcDzwtrbrz8fcMt79m1mh9Z9Y/RbAiAICZEc5R9sGOw3ryL7sMt+fZMgAk\nH8I5Sqq/qNPqF7cabj94ULpW3H4+85YBIAnxyR9hfb2FLUllC6dp3LCcCFUEAIg3hHME9TWYbyrJ\n10XnjIxgRQCAeEQ4R8g/jjSI17A9AAAK8klEQVTokT98ZKitPV1acydLbwIATiEN+uhIbZN+8ttN\nChicI/XIrTM0YkhWZIsCAMQ1wtmg1rY2lT+9WcfqWgy1t1qk3/zvi1hMBADwDYSzAX19tjxsSKbK\nF00nmAEAXSKcw9Tc6jMczBPOzNa//HOBcuwDIlwVACCREM4hamxq1Ue7jum5N/YYar9m8UzlZWVE\nuCoAQCIinHvR2tamR579SIfcJw21n3HWEN1+1WRuYQMAQkY498Dr82v505t1rN7YoK/HfliknIEs\nvQkACA/h3AV/IKAX3titjZ8ekc8ffvsMm1VrFl+gjPS0yBcHAEh4hPPX+AMBPfTMhzroCv82dlaG\nVctumiFnbmYUKgMAJAvC+WteeKPaUDBfOPUM3XTZWTxbBgD0GeH8pSavT8+vr9YHnx0Lq12KpJ/d\ndaGyMriFDQCIjKQP5yZvm15+s1pbdtfI6wuE1daZm66Hb2FbRwBAZCVtqvgDAVW8vVfvfXJYLa2h\nh7I1RZowKke3XlWgHDsjsQEAkZeU4ez1+fX8+t3auONoWO2G5mXoxwunK9PGLWwAQPQkVTj7AwG9\n9Ga1tu5xq97TGnK79FSLLvj2UF1fMoEBXwCAqEuacPb7A3r42S36osYTVrvzJzm18IqzZEuzRqky\nAAA6MxTOPp9PS5cu1eHDh2W1WrVq1SqNGjWq03tee+01/f73v1dKSopmzpype+65JyIFh8vr88tV\n16RnXt8VVjAPSLdq1pRhmjd3HFfLAIB+ZSic//KXv2jQoEFau3at3nvvPa1du1aPP/54x/Hm5mat\nWbNG69at08CBA3Xdddfpqquu0rhx4yJWeG/8gYBe/usebfz0SFgDvmxpKZo+wan5JfnKtCXNjQUA\ngIkYSp+qqipdffXVkqQLLrhAZWVlnY5nZGRo3bp1stvtkqScnBzV19f3sdTwVLy9V29/dCjk9+fa\n03XW6DwtKBnPgC8AQEwZCme32628vDxJUkpKiiwWi1pbW5Went7xnvZg3r17tw4dOqSpU6dGoNzQ\neH1+fby7JuT3nz/pDC28YiLPlQEAptBrOL/66qt69dVXO31v+/btnV4Hg8Eu2+7fv19LlizR2rVr\nlZbW89Vobm6mUlMjE45H3CdV2xjaaOx/Gj5ISxedK6uV58qncziyYl1CXKCfQkM/hYZ+Ck0y9FOv\n4Xzttdfq2muv7fS9pUuXyuVyaeLEifL5fAoGg52umiXp6NGjWrx4sR599FGdddZZvRZSV9cUZund\n8/v8ystK7zGgcwamq3CCQwuKx6u21thezYnK4ciSy9UY6zJMj34KDf0UGvopNInUTz39kWHocrGo\nqEivv/66JGnDhg0677zzvvGeZcuWafny5SooKDDyI/rElmbVtAnObo8XTR6qVXfM1I2XMm8ZAGA+\nhp45X3nlldq4caPmz5+v9PR0rV69WpL01FNPacaMGcrJydGWLVv085//vKPNokWLdMkll0Sm6hDM\nmztOgWBQGz89qpbWU5syD0i3qujbQ1V6yXhCGQBgWpZgdw+M+1m0blO0z3POzRuo1GCQQV8hSKTb\nRtFEP4WGfgoN/RSaROqniN/Wjie2NKtGOrM0elg2wQwAiAsJH84AAMQbwhkAAJMhnAEAMBnCGQAA\nkyGcAQAwGcIZAACTIZwBADAZwhkAAJMhnAEAMBnCGQAAkyGcAQAwGdNsfAEAAE7hyhkAAJMhnAEA\nMBnCGQAAkyGcAQAwGcIZAACTIZwBADCZ1FgXEC0+n09Lly7V4cOHZbVatWrVKo0aNarTe1577TX9\n/ve/V0pKimbOnKl77rknRtXGxsqVK7V9+3ZZLBaVlZVpypQpHcc2btyoxx57TFarVbNnz9bixYtj\nWGls9dRPmzZt0mOPPaaUlBSNGTNGK1asUEpKcv7N21M/tVu7dq22bdum559/PgYVmkNP/XTkyBHd\ne++98vl8mjRpkh5++OEYVhpbPfXTiy++qHXr1iklJUWTJ0/WsmXLYlhplAQTVGVlZXD58uXBYDAY\nfPfdd4N33XVXp+NNTU3Biy++ONjY2BgMBALB733ve8E9e/bEotSY+OCDD4K33357MBgMBvfu3Ru8\n7rrrOh2/4oorgocPHw76/f7g/Pnzk6pvTtdbP5WUlASPHDkSDAaDwR/+8IfBd955p99rNIPe+ikY\nDAb37NkTnDdvXvCGG27o7/JMo7d+uvPOO4NvvPFGMBgMBpcvXx48dOhQv9doBj31U2NjY/Diiy8O\n+ny+YDAYDN58883BrVu3xqTOaErYP/GrqqpUUlIiSbrgggv08ccfdzqekZGhdevWyW63y2KxKCcn\nR/X19bEoNSaqqqpUXFwsSRo7dqwaGhrk8XgkSV988YWys7M1bNgwpaSkaM6cOaqqqopluTHTUz9J\nUmVlpYYOHSpJysvLU11dXUzqjLXe+kmSVq9enXR3p76up34KBAL66KOPNHfuXElSeXm5hg8fHrNa\nY6mnfkpLS1NaWpqamprU1tam5uZmZWdnx7LcqEjYcHa73crLy5MkpaSkyGKxqLW1tdN77Ha7JGn3\n7t06dOiQpk6d2u91xorb7VZubm7H67y8PLlcLkmSy+Xq6LuvH0s2PfWT9NXvUE1Njd5//33NmTOn\n32s0g976qbKyUueee65GjBgRi/JMo6d+qq2t1cCBA7Vq1SrNnz9fa9eujVWZMddTP9lsNi1evFjF\nxcW6+OKLNXXqVI0ZMyZWpUZNQjxzfvXVV/Xqq692+t727ds7vQ52s0rp/v37tWTJEq1du1ZpaWlR\nq9HsuusfdNZVPx0/flx33HGHysvLO32gJLPT+6m+vl6VlZV65plndOzYsRhWZT6n91MwGNSxY8d0\n0003acSIEbr99tv1zjvv6KKLLopdgSZxej95PB49+eSTev3112W327Vw4ULt2rVLEydOjGGFkZcQ\n4Xzttdfq2muv7fS9pUuXyuVyaeLEifL5fAoGg0pPT+/0nqNHj2rx4sV69NFHddZZZ/VnyTHndDrl\ndrs7XtfU1MjhcHR57NixY3I6nf1eoxn01E/SqQ+K2267TXfffbdmzZoVixJNoad+2rRpk2pra3X9\n9dertbVVn3/+uVauXKmysrJYlRszPfVTbm6uhg8frjPPPFOSNHPmTO3Zsycpw7mnftq3b59GjRrV\ncXdv+vTp2rFjR8KFc8Le1i4qKtLrr78uSdqwYYPOO++8b7xn2bJlWr58uQoKCvq7vJgrKirS+vXr\nJUk7d+6U0+nsuEU7cuRIeTweHTx4UG1tbdqwYYOKiopiWW7M9NRP0qnnqAsXLtTs2bNjVaIp9NRP\nl19+uV577TX98Y9/1C9+8QsVFBQkZTBLPfdTamqqRo0apf3793ccT8TbtaHoqZ9GjBihffv2qaWl\nRZK0Y8cOjR49OlalRk3C7krl9/v14x//WPv371d6erpWr16tYcOG6amnntKMGTOUk5Ojq6++utPw\n/EWLFumSSy6JYdX9a82aNdqyZYssFovKy8v12WefKSsrSyUlJfrwww+1Zs0aSdKll16qW265JcbV\nxk53/TRr1izNmDFDhYWFHe/97ne/q3nz5sWw2tjp6fep3cGDB/XAAw8k9VSqnvrpwIEDWrp0qYLB\noPLz87V8+fKknZrXUz+98sorqqyslNVqVWFhoe6///5YlxtxCRvOAADEq+T8kwwAABMjnAEAMBnC\nGQAAkyGcAQAwGcIZAACTIZwBADAZwhkAAJMhnAEAMJn/D3v9Is9HUWMDAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 576x396 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "metadata": {
        "id": "-dw9TGERe3cw",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## Useful input features for NN (Neural Networks)\n",
        "\n",
        "The \"apc\" values are typically used as input to the NN for contact cleaning or structure prediction. Though in recent advances (aka DeepMind/Alphafold), the entire MRF was used as the input. More specificially LxLx442. The 442 channels are the 21x21 + (raw and/or apc) value."
      ]
    },
    {
      "metadata": {
        "id": "nqbDMKdsYG7Q",
        "colab_type": "code",
        "outputId": "60e40a12-b1ad-4ea4-b8bd-5fa757ff8a48",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 51
        }
      },
      "cell_type": "code",
      "source": [
        "w_out = np.zeros((msa[\"ncol_ori\"], msa[\"ncol_ori\"], states * states + 1))\n",
        "v_out = np.zeros((msa[\"ncol_ori\"], states))\n",
        "\n",
        "mrf_ = np.reshape(mrf[\"w\"],(-1, states * states))\n",
        "mtx_ = np.expand_dims(mtx[\"apc\"],-1)\n",
        "\n",
        "w_out[(mtx[\"i\"],mtx[\"j\"])] = np.concatenate((mrf_, mtx_),-1)\n",
        "w_out += np.transpose(w_out,(1,0,2))\n",
        "v_out[mrf[\"v_idx\"]] = mrf[\"v\"]\n",
        "\n",
        "print(\"w_out\",w_out.shape)\n",
        "print(\"v_out\",v_out.shape)"
      ],
      "execution_count": 17,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "w_out (62, 62, 442)\n",
            "v_out (62, 21)\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}